RHESSysIOinR_Coweeta_WS32_74
Carlos Quintero
6/21/2023
Coweeta Hydrologic Laboratory Watershed 32 using Static, SSURGO, and DSM soil inputs. Previously WS27.
The following R Markdown script prepares and runs the Regional Hydro-Ecologic Simulation System (RHESSYS) within the R environment. Modified from RHESSysIOinR example scripts created by Will Burke and Ryan Bart. Includes CLHS code snippets and SDA code from Dylan Beaudette.
Windows Subsystem via Linux must be installed before RHESSys will run in a Windows environment. GRASS8 must be installed to run the spatial preprocessing portion of the script.
Code has been modified to run RHESSys 7.4 and 5.18. 7.4 vegetation processing differs slightly from 5.18. Various template differences exist between the two versions of RHESSys and template files must be carefully reviewed if moving between versions. This script is currently set up to run RHESSys 7.4
#install.packages("devtools")
#devtools::install_github("ropensci/FedData")
library(aqp)## This is aqp 1.41library(daymetr)
library(clhs)## Registered S3 methods overwritten by 'tibble':
## method from
## format.tbl pillar
## print.tbl pillarlibrary(EcoHydRology)## Loading required package: operators##
## Attaching package: 'operators'## The following objects are masked from 'package:base':
##
## options, strrep## Loading required package: topmodel## Loading required package: DEoptim## Loading required package: parallel##
## DEoptim package
## Differential Evolution algorithm in R
## Authors: D. Ardia, K. Mullen, B. Peterson and J. Ulrich## Loading required package: XMLlibrary(elevatr)
library(ezknitr)
library(FedData)##
## Attaching package: 'FedData'## The following object is masked from 'package:operators':
##
## %>%library(foreach)
library(foreign)
library(ggplot2)
library(hydroGOF)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Attaching package: 'hydroGOF'## The following object is masked from 'package:topmodel':
##
## NSefflibrary(imputeTS)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Attaching package: 'imputeTS'## The following object is masked from 'package:zoo':
##
## na.locf## The following object is masked from 'package:operators':
##
## %>%library(knitr)
library(latticeExtra)## Loading required package: lattice##
## Attaching package: 'latticeExtra'## The following object is masked from 'package:ggplot2':
##
## layerlibrary(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(Metrics)##
## Attaching package: 'Metrics'## The following objects are masked from 'package:hydroGOF':
##
## mae, mse, rmselibrary(raster)## Loading required package: sp##
## Attaching package: 'raster'## The following objects are masked from 'package:aqp':
##
## metadata, metadata<-library(rasterVis)
library(readxl)
library(reshape2)
library(rgdal)## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.1, released 2021/12/27
## Path to GDAL shared files: C:/Users/Carlos/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/Carlos/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-6
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.library(rgrass)## GRASS GIS interface loaded with GRASS version: (GRASS not running)library(Rmisc)## Loading required package: plyrlibrary(RHESSysPreprocessing)
library(RHESSysIOinR)##
## Attaching package: 'RHESSysIOinR'## The following object is masked from 'package:RHESSysPreprocessing':
##
## read_worldlibrary(sf)## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE##
## Attaching package: 'sf'## The following object is masked from 'package:operators':
##
## %>%library(soilDB)
library(sp)
library(stringi)
library(tactile)
library(terra)## terra 1.5.21##
## Attaching package: 'terra'## The following object is masked from 'package:rgdal':
##
## project## The following object is masked from 'package:knitr':
##
## spin## The following object is masked from 'package:zoo':
##
## time<-## The following object is masked from 'package:ggplot2':
##
## arrowlibrary(testthat)##
## Attaching package: 'testthat'## The following object is masked from 'package:terra':
##
## describe## The following object is masked from 'package:operators':
##
## %>%library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v tibble 3.0.3 v dplyr 1.0.5
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1
## v purrr 0.3.4## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x forcats::%>%() masks stringr::%>%(), dplyr::%>%(), purrr::%>%(), tidyr::%>%(), tibble::%>%(), testthat::%>%(), sf::%>%(), imputeTS::%>%(), FedData::%>%(), operators::%>%()
## x purrr::accumulate() masks foreach::accumulate()
## x dplyr::arrange() masks plyr::arrange()
## x terra::arrow() masks ggplot2::arrow()
## x lubridate::as.difftime() masks base::as.difftime()
## x dplyr::combine() masks aqp::combine()
## x purrr::compact() masks plyr::compact()
## x dplyr::count() masks plyr::count()
## x lubridate::date() masks base::date()
## x readr::edition_get() masks testthat::edition_get()
## x tidyr::extract() masks terra::extract(), raster::extract()
## x dplyr::failwith() masks plyr::failwith()
## x dplyr::filter() masks stats::filter()
## x dplyr::id() masks plyr::id()
## x terra::intersect() masks raster::intersect(), lubridate::intersect(), base::intersect()
## x purrr::is_null() masks testthat::is_null()
## x dplyr::lag() masks stats::lag()
## x latticeExtra::layer() masks ggplot2::layer()
## x readr::local_edition() masks testthat::local_edition()
## x dplyr::matches() masks tidyr::matches(), testthat::matches()
## x dplyr::mutate() masks plyr::mutate()
## x dplyr::rename() masks plyr::rename()
## x dplyr::select() masks raster::select()
## x lubridate::setdiff() masks base::setdiff()
## x dplyr::slice() masks aqp::slice()
## x dplyr::src() masks terra::src()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
## x terra::union() masks raster::union(), lubridate::union(), base::union()
## x purrr::when() masks foreach::when()library(viridisLite)
knitr::opts_knit$set(root.dir = system.file("extdata/", package = "RHESSysIOinR"))
## Prep RHESSys
expect_file_exists = function(path) {
# 1. Capture object and label
act <- quasi_label(rlang::enquo(path), arg = "path")
# 2. Call expect()
expect(
file.exists(act$val),
sprintf("%s file does not exist.", act$lab)
)
# 3. Invisibly return the value
invisible(act$val)
}
expect_file_sizeKB_gt = function(path, size_KB) {
# 1. Capture object and label
act <- quasi_label(rlang::enquo(path), arg = "path")
# 2. Call expect()
expect(
(file.size(act$val) / 1024) > size_KB,
sprintf("%s file is not greater than %f KB.", act$lab, size_KB)
)
# 3. Invisibly return the value
invisible(act$val)
}
setwd(system.file("extdata/", package = "RHESSysIOinR"))
### Downloads most recent version of RHESSYs
# gert::git_clone(url = "https://github.com/RHESSys/RHESSys", path = "./rh_dev", branch = "master")
### Make sure to change files in rh_dev to 5.18 if using older RHESSys files
# compile_rhessys(location = "rh_dev/")
rh_ver = dir(path = "rh_dev/rhessys/", pattern = "^rhessys\\d+",recursive = F)
test_that("compile_rhessys works + rhessys compiles via R system", {
expect_gt(length(rh_ver), 0)
})## Test passedrh_path = file.path("rh_dev/rhessys/", rh_ver)Load in Shapefiles for Watershed Setup
getwd()## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"setwd(system.file("extdata/", package = "RHESSysIOinR"))
watershedshapefile <-("C:/Users/Carlos/Desktop/RHESSys_R_Processing/Boundary Shapefiles/Coweeta_Hydrologic_Laboratory.shp")
streamshapefile<- ("C:/Users/Carlos/Desktop/ORISE/GIS/CW/Streamflow/StreamsCaldwellClipped.shp")
weirshapefile<- ("C:/Users/Carlos/Desktop/ORISE/GIS/CW/Weirs/Subbasinoutlets.shp")
roadshapefile <- ( "C:/Users/Carlos/Desktop/ORISE/GIS/CW/Roads/UpdatedRoads4Reprojected32617.shp")
#streamflow dataset
#soil moisture dataset##Quick way to create base file necessary for climate processing not functional yetDownload Daymet Data for Watershed of Interest using daymetr package Currently only creates precip, tmin and tmax files.
Extend DAYMET datasets to allow for longer spinup period Read Daymet Datasets and Replicate 5x times to create a longer period of record for spin up Create base File manually for now
## this may be the source of lag in outputs/measured data
getwd()## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"precip <- read.csv('clim/cwtws32daymet.rain', header = FALSE)
tmin <- read.csv('clim/cwtws32daymet.tmin', header = FALSE)
tmax <- read.csv('clim/cwtws32daymet.tmax', header = FALSE)
#tavg <- read.csv('clim/cwtws32daymet.tavg', header = FALSE)
#rh <- read.csv('clim/cwtws32daymet.relative_humidity', header = FALSE)
#wind <- read.csv('clim/cwtws32daymet.wind', header = FALSE)
head(precip)## V1
## 1 1980 1 1 1
## 2 0.00427
## 3 0
## 4 0.00866
## 5 0.01468
## 6 0nrow(precip)-1## [1] 15341substr(precip[1,1],1,nchar(precip[1,1])-2)## [1] "1980 1 1"dt <- as.Date(substr(precip[1,1],1,nchar(precip[1,1])-2), "%Y %m %d")
dt.new<- dt - as.difftime((nrow(precip)-1)*5, units="days")
#remove first value and repeats dataset 6 times
longcwtraindata<- do.call("rbind",replicate(6,as.list((precip[-1,]))))
longcwttmindata<- do.call("rbind",replicate(6,as.list((tmin[-1,]))))
longcwttmaxdata<- do.call("rbind",replicate(6,as.list((tmax[-1,]))))
#longcwttavgdata<- do.call("rbind",replicate(6,as.list((tavg[-1,]))))
#longcwtrhdata<- do.call("rbind",replicate(6,as.list((rh[-1,]))))
#longcwtwinddata<- do.call("rbind",replicate(6,as.list((wind[-1,]))))
longcwstartdate<- paste(year(ymd(dt.new)),month(ymd(dt.new)),day(ymd(dt.new)),"1")
longcwtrain<- rbind(longcwstartdate,longcwtraindata)
longcwttmin<- rbind(longcwstartdate,longcwttmindata)
longcwttmax<- rbind(longcwstartdate,longcwttmaxdata)
#longcwttavg<- rbind(longcwstartdate,longcwttavgdata)
#longcwtrh<- rbind(longcwstartdate,longcwtrhdata)
#longcwtwind<- rbind(longcwstartdate,longcwtwinddata)
longcwtrain<- data.frame(longcwtrain)
longcwttmin<- data.frame(longcwttmin)
longcwtmax<- data.frame(longcwttmax)
#longcwttavg<- data.frame(longcwttavg)
#longcwtrh<- data.frame(longcwtrh)
#longcwtwind<- data.frame(longcwtwind)
write.table(longcwtrain, file = "clim/cwtws32extended.rain", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(longcwttmin, file = "clim/cwtws32extended.tmin", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(longcwttmax, file = "clim/cwtws32extended.tmax", quote = FALSE, row.names = FALSE, col.names = FALSE)
#write.table(longcwttavg, file = "clim/cwtws32extended.tavg", quote = FALSE, row.names = FALSE, col.names = FALSE)
#write.table(longcwtrh, file = "clim/cwtws32extended.relative_humidity", quote = FALSE, row.names = FALSE, col.names = FALSE)
#write.table(longcwtwind, file = "clim/cwtws32extended.wind", quote = FALSE, row.names = FALSE, col.names = FALSE)
read_clim("clim/cwtws32extended.rain", dates_out=TRUE)## [1] "1769 12 27 01" "2021 12 31 24"read_clim("clim/cwtws32daymet.rain", dates_out=TRUE)## [1] "1980 01 01 01" "2021 12 31 24"#daymet_data$data
# precip data for watershed not readily availableDownload and Prepare DEM Code to download DEM originally prepared by Dylan Beaudette
# use watershed polygons for BBOX to request elevation data
x <- read_sf(watershedshapefile)
# buffer to ensure that there are no truncated watersheds
x.buff <- st_as_sf(st_buffer(st_as_sfc(st_bbox(x)), dist = 1000))
# convert to GCS WGS84 -> no transformation / warp will be applied to DEM
x.gcs <- st_transform(x, 4326)
x.buff.gcs <- st_transform(x.buff, 4326)
# requires sf collection (geometry + attributes)
# get DEM in GCS WGS84 and use z = 14 for best available data
e <- get_elev_raster(locations = x.buff.gcs, z = 14, clip = 'bbox')## Mosaicing & Projecting## Clipping DEM to bbox## Note: Elevation units are in meters.res(e)## [1] 3.921794e-05 3.921794e-05# check: OK
{plot(e, main = "DEM, Watershed Outline and Buffer" )
plot(st_geometry(x.buff.gcs), add = TRUE, col = NA, lwd = 2)
plot(st_geometry(x.gcs), add = TRUE, col = NA)}
{plot(e, main = "DEM and selected Watershed")
plot(st_geometry(x.gcs[which(x$WS=='32'),]), add = TRUE, col = NA)}
cropped_dem<- crop(e,st_bbox(x.gcs))
plot(cropped_dem, main = "DEM cropped to Buffer")
plot(st_geometry(x.gcs[which(x$WS=='32'),]), add = TRUE, col = NA)
newproj<- "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"
reproje<- projectRaster(from = cropped_dem, crs = newproj, method = 'ngb')
plot(reproje, main = "Reprojected cropped DEM")
dummydata<- reproje
res(dummydata) <- 10
resampleddem <- resample(reproje, dummydata, method = 'bilinear')
plot(resampleddem, main = "DEM reprojected to 10m Resolution")
Start Spatial Processing in GRASS8 Code to prepare spatial input files for RHESSYs takes cues from Will Burkes “spatial_input_gen.R” script but is prepared for GRASS8 instead of GRASS7
#gisbase_grass <- ifelse(.Platform$OS.type == "windows", link2GI::paramGRASSw()$gisbase_GRASS[1], link2GI::paramGRASSx()$gisbase_GRASS[1])
## set easting and northing as coordinates of lowest point within boundary basin, set threshold for basin delineation
threshold = 400
easting = NULL
northing = NULL
initGRASS(gisBase = "C:/Program Files/GRASS GIS 8.2",
gisDbase = "C:/Users/Carlos/Documents/grassdata",
location = "Coweeta",
mapset = "PERMANENT",
override = TRUE)## Warning in system(syscmd, intern = intern, ignore.stderr = ignore.stderr, :
## running command 'g.proj.exe -w' had status 884## gisdbase C:/Users/Carlos/Documents/grassdata
## location Coweeta
## mapset PERMANENT
## rows 611
## columns 537
## north 3884060
## south 3877950
## west 273846.7
## east 279216.7
## nsres 10
## ewres 10
## projection +proj=utm +no_defs +zone=17 +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000 +type=crs +to_meter=1execGRASS("g.proj", proj4="+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs", flag = "c" )
execGRASS("g.gisenv", flag="n")
execGRASS("g.region", flag=c("p"))
##Convert DEM to Spatrast and then import DEM to GRASS
spatrast_dem10 <- rast(resampleddem)
write_RAST(spatrast_dem10, "grass_dem10", overwrite = TRUE)## SpatRaster read into GRASS using r.in.gdal from memory## The following code is heavily influenced by Will Burkes "spatial_input_gen.R" script
execGRASS("g.region", raster = "grass_dem10", flag=c("p"))
## Import shapefiles (Streams, Weirs, Watersheds, Roads)
## only roads and Watersheds necessary to Run
streams<- vect("C:/Users/Carlos/Desktop/ORISE/GIS/CW/Streamflow/StreamsCaldwellClipped.shp")
streams<- project(streams,"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs")
write_VECT(streams, "streams")
weirs<- vect("C:/Users/Carlos/Desktop/ORISE/GIS/CW/Weirs/Subbasinoutlets.shp")
weirs<- project(weirs,"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs")
write_VECT(weirs, "weirs")
watersheds<- vect("C:/Users/Carlos/Desktop/RHESSys_R_Processing/Boundary Shapefiles/Coweeta_Hydrologic_Laboratory.shp")
wgs84watersheds<- project(watersheds,"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs")
write_VECT(wgs84watersheds, "watersheds")
roads <- vect( "C:/Users/Carlos/Desktop/ORISE/GIS/CW/Roads/UpdatedRoads4Reprojected32617.shp")## Warning: [vect] Z coordinates ignoredroads<- project(roads,"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs")
write_VECT(roads,"roads")
## Check all vectors that have been exported to GRASS
execGRASS("g.list", parameters = list(type = "vector"))
## Generate a depressionless DEM and d-8 flow direction maps from DEM
execGRASS("r.fill.dir", input="grass_dem10", output="dem10f", direction="dir10", flags = c("overwrite") )
## Generate flow accumulation, drainage, and horizon maps
execGRASS("r.watershed", elevation = "dem10f", accumulation = "acc10", drainage = "drain10")
## Generate slope, aspect, wetness index, horizon, and road maps
execGRASS("r.slope.aspect", elevation = "grass_dem10", slope = "slope10", aspect = "aspect10")
execGRASS("r.topidx", input = "grass_dem10", output = "topidx10")
execGRASS("r.mapcalc", expression = "topidx10_100 = topidx10 * 100", flags = "overwrite")
execGRASS("r.horizon", flags = "d", elevation = "grass_dem10", direction = 0, output = "east")
execGRASS("r.horizon", flags = "d", elevation = "grass_dem10", direction = 180, output = "west")
execGRASS("r.mapcalc", expression = "e_horizon_1000 = sin(east_000) * 1000", flags = "overwrite")
execGRASS("r.mapcalc", expression = "w_horizon_1000 = sin(west_180) * 1000", flags = "overwrite")
## Generate a hillslope and stream map
subbasin_name = paste0("subbasin_th",as.character(threshold))
stream_name = paste0("stream_th",as.character(threshold))
hillslope_name = paste0("hillslope_th",as.character(threshold))
execGRASS("r.watershed", elevation = "dem10f", threshold = threshold, basin = subbasin_name, stream = stream_name, half_basin = hillslope_name)
## Delineate watershed boundary
weirs[12,]## class : SpatVector
## geometry : points
## dimensions : 1, 1 (geometries, attributes)
## extent : 276120.4, 276120.4, 3881339, 3881339 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=17 +datum=WGS84 +units=m +no_defs
## names : id
## type : <int>
## values : 12easting = 276122.49
northing = 3881345.46
execGRASS("g.region", raster = "grass_dem10", flag=c("p"))
# r.water.outlet input=drain10 output=basin_ws321 coordinates=276122.49,3881345.46 --overwrite
execGRASS("r.water.outlet", input ="drain10", output="basin_ws32",coordinates = c(easting, northing), flags = "overwrite")
## Mask Data
# execGRASS("r.mask", input = "basin_ws32", maskcats = 1)
execGRASS("r.info",map = "basin_ws32")
execGRASS("r.info",map = "grass_dem10")
## Generate Patch Maps
execGRASS("r.mapcalc", expression = "patch = row() * 537 + col()", flags = "overwrite")
## Create Road Maps
execGRASS("v.to.rast", input="roads" ,output="roads", use="cat", label_column="ROADNAME", flags = "overwrite")
execGRASS("r.mapcalc", expression = "roads = roads > 0", flags = "overwrite")
execGRASS("r.null", map = "roads", null= 0)
## Check all rasters that have been exported to or created in GRASS
execGRASS("g.list", parameters = list(type = "raster"))
## Export Prepared Maps from GRASS to R
filleddem <- read_RAST("dem10f")
##Export maps back into R or export directly from GRASS to disk
grass_acc10 <- read_RAST("acc10")
grass_aspect10 <- read_RAST("aspect10")
grass_basin_ws32 <- read_RAST("basin_ws32")
grass_dem10 <- read_RAST("grass_dem10")
grass_dem10f <- read_RAST("dem10f")
grass_dir10 <- read_RAST("dir10")
grass_drain10 <- read_RAST("drain10")
grass_e_horizon_1000 <- read_RAST("e_horizon_1000")
grass_hillslope_th1000 <- read_RAST("hillslope_th1000")
#grass_patch <- read_RAST("patch")
grass_slope10 <- read_RAST("slope10")
grass_stream_th1000 <- read_RAST("stream_th1000")
grass_topidx10 <- read_RAST("topidx10")
grass_topidx10_100 <- read_RAST("topidx10_100")
grass_w_horizon_1000 <- read_RAST("w_horizon_1000")
grass_east_000 <- read_RAST("east_000")
grass_west_180 <- read_RAST("west_180")
grass_subbasin_th1000 <- read_RAST("subbasin_th1000")
execGRASS("r.out.gdal", input = "acc10", output = "grassexport/cwt/acc10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "aspect10", output = "grassexport/cwt/aspect10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "basin_ws32", output = "grassexport/cwt/basin_ws32.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "grass_dem10", output = "grassexport/cwt/grass_dem10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "dem10f", output = "grassexport/cwt/dem10f.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "dir10", output = "grassexport/cwt/dir10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "drain10", output = "grassexport/cwt/drain10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "e_horizon_1000", output = "grassexport/cwt/e_horizon_1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "hillslope_th1000", output = "grassexport/cwt/hillslope_th1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "patch", output = "grassexport/cwt/patch.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "slope10", output = "grassexport/cwt/slope10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "stream_th1000", output = "grassexport/cwt/stream_th1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "topidx10", output = "grassexport/cwt/topidx10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "topidx10_100", output = "grassexport/cwt/topidx10_100.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "w_horizon_1000", output = "grassexport/cwt/w_horizon_1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "east_000", output = "grassexport/cwt/east_000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "west_180", output = "grassexport/cwt/west_180.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "subbasin_th1000", output = "grassexport/cwt/subbasin_th1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "roads", output = "grassexport/cwt/roads.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "stream_th1500", output = "grassexport/cwt/stream_th1500.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "stream_th500", output = "grassexport/cwt/stream_th500.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "stream_th400", output = "grassexport/cwt/stream_th400.tif", format = "GTiff", flags = "overwrite")Crop Maps that were prepared in GRASS. Maps are cropped R because MASK command is not functioning properly for GRASS via R. Original input maps are been cropped before processing worldfile and flowtable in R, without cropping the RHESSYs worldfile and flowtable preparation in the R environment will not function properly,
ws32<- raster("grassexport/cwt/basin_ws32.tif")
acc10<- raster("grassexport/cwt/acc10.tif")
patches<- raster("grassexport/cwt/patch.tif")
aspect10<- raster("grassexport/cwt/aspect10.tif")
dem10<- raster("grassexport/cwt/grass_dem10.tif")
dem10f<- raster("grassexport/cwt/dem10f.tif")
dir10<- raster("grassexport/cwt/dir10.tif")
drain10<- raster("grassexport/cwt/drain10.tif")
hillslope<- raster("grassexport/cwt/hillslope_th1000.tif")
roads<- brick("grassexport/cwt/roads.tif")
slope10<- raster("grassexport/cwt/slope10.tif")
streams<- raster("grassexport/cwt/stream_th1000.tif")
topidx10_100<- raster("grassexport/cwt/topidx10_100.tif")
e_horizon_1000 <- raster("grassexport/cwt/e_horizon_1000.tif")
w_horizon_1000 <- raster("grassexport/cwt/w_horizon_1000.tif")
subbasins<- raster("grassexport/cwt/subbasin_th1000.tif")
ws32polygon<- rasterToPolygons(ws32, fun = function(x){x==1}, dissolve = TRUE, digits = 12)## Loading required namespace: rgeosmaskedws32 <- crop(mask(ws32,ws32polygon),ws32polygon)
maskedacc10 <- crop(mask(acc10,ws32polygon),ws32polygon)
maskedpatches <- crop(mask(patches,ws32polygon),ws32polygon)
maskedaspect10 <- crop(mask(aspect10,ws32polygon),ws32polygon)
maskeddem10 <- crop(mask(dem10,ws32polygon),ws32polygon)
maskeddem10f <- crop(mask(dem10f,ws32polygon),ws32polygon)
maskeddir10 <- crop(mask(dir10,ws32polygon),ws32polygon)
maskeddrain10 <- crop(mask(drain10,ws32polygon),ws32polygon)
maskedhillslope <- crop(mask(hillslope,ws32polygon),ws32polygon)
maskedroads <- crop(mask(roads,ws32polygon),ws32polygon)
maskedslope10 <- crop(mask(slope10,ws32polygon),ws32polygon)
maskedstreams <- crop(mask(streams,ws32polygon),ws32polygon)
maskedtopidx10_100 <- crop(mask(topidx10_100,ws32polygon),ws32polygon)
maskede_horizon_1000 <- crop(mask(e_horizon_1000,ws32polygon),ws32polygon)
maskedw_horizon_1000 <- crop(mask(w_horizon_1000,ws32polygon),ws32polygon)
maskedsubbasins<- crop(mask(subbasins,ws32polygon),ws32polygon)
## Save as tif
writeRaster(maskedws32,"grassexport/cwt_ws32/basin_ws32",format = "GTiff", overwrite = TRUE)
writeRaster(maskedacc10,"grassexport/cwt_ws32/acc10",format = "GTiff", overwrite = TRUE)
writeRaster(maskedpatches,"grassexport/cwt_ws32/patch",format = "GTiff", overwrite = TRUE)
writeRaster(maskedaspect10,"grassexport/cwt_ws32/aspect10",format = "GTiff", overwrite = TRUE)
writeRaster(maskeddem10,"grassexport/cwt_ws32/dem10",format = "GTiff", overwrite = TRUE)
writeRaster(maskeddem10f,"grassexport/cwt_ws32/dem10f",format = "GTiff", overwrite = TRUE)
writeRaster(maskeddir10,"grassexport/cwt_ws32/dir10",format = "GTiff", overwrite = TRUE)
writeRaster(maskeddrain10,"grassexport/cwt_ws32/drain10",format = "GTiff", overwrite = TRUE)
writeRaster(maskedhillslope,"grassexport/cwt_ws32/hillslope",format = "GTiff", overwrite = TRUE)
writeRaster(maskedroads,"grassexport/cwt_ws32/roads",format = "GTiff", overwrite = TRUE)
writeRaster(maskedslope10,"grassexport/cwt_ws32/slope10",format = "GTiff", overwrite = TRUE)
writeRaster(maskedstreams,"grassexport/cwt_ws32/streams",format = "GTiff", overwrite = TRUE)
writeRaster(maskedtopidx10_100,"grassexport/cwt_ws32/topidx10_100",format = "GTiff", overwrite = TRUE)
writeRaster(maskede_horizon_1000,"grassexport/cwt_ws32/e_horizon_1000",format = "GTiff", overwrite = TRUE)
writeRaster(maskedw_horizon_1000,"grassexport/cwt_ws32/w_horizon_1000",format = "GTiff", overwrite = TRUE)
writeRaster(maskedsubbasins,"grassexport/cwt_ws32/subbasins",format = "GTiff", overwrite = TRUE)
plot(maskedws32, main ="ws32")
plot(maskedpatches, main = "patches")
plot(maskedacc10, main = "acc10")
plot(maskedaspect10, main = "aspect10")
plot(maskedsubbasins, main = "subbasins")
plot(maskeddem10, main = "dem10")
plot(maskeddem10f, main = "dem10f")
plot(maskeddir10, main = "dir10")
plot(maskeddrain10, main = "drain10")
plot(maskedhillslope, main = "hillslope")
plot(maskedroads, main = "roads")
plot(maskedslope10, main = "slope10")
plot(maskedstreams, main = "streams")
plot(maskedtopidx10_100, main = "topidx10_100")
plot(maskede_horizon_1000, main = "e_horizon_1000")
plot(maskedw_horizon_1000, main = "w_horizon_1000")
SSURGO-via-SDA Code originally prepared by Dylan Beaudette Download SSURGO data and make sure CRS/Extent/Crop are Correct for input into RHESSys
# load BBOX
x <- read_sf(watershedshapefile)
# get gSSURGO grid here
mu <- mukey.wcs(aoi = x, db = 'gssurgo')
# unique map unit keys
ll <- levels(mu)[[1]]
# note SSA boundary
levelplot(mu, att = 'ID', margin = FALSE, colorkey = FALSE, col.regions = viridis, main = "SSURGO Map Units")
## Note: some of these map units are dominated by non-soil components
# Dominant Component (Numeric) -> NODATA
# Weighted Average -> use any available data, but wt. mean are diluted (see source data)
# get thematic data from SDA
# dominant component
# depth-weighted average
# sand, silt, clay, AWC (RV)
p <- get_SDA_property(property = c("sandtotal_r","silttotal_r","claytotal_r", "awc_r"),
method = "Weighted Average",
mukeys = ll$ID,
top_depth = 0,
bottom_depth = 200)
head(p)## areasymbol musym
## 1 NC113 BuD
## 2 NC113 BuF
## 3 NC113 CaF
## 4 NC113 CdD
## 5 NC113 CdE
## 6 NC113 CdF
## muname
## 1 Burton-Craggey-Rock outcrop complex, windswept, 15 to 30 percent slopes, stony
## 2 Burton-Craggey-Rock outcrop complex, windswept, 30 to 95 percent slopes, stony
## 3 Cashiers gravelly fine sandy loam, 50 to 95 percent slopes
## 4 Chandler gravelly fine sandy loam, 15 to 30 percent slopes
## 5 Chandler gravelly fine sandy loam, 30 to 50 percent slopes
## 6 Chandler gravelly fine sandy loam, 50 to 95 percent slopes
## mukey sandtotal_r silttotal_r claytotal_r awc_r
## 1 545800 66.36 19.50 14.14 0.05
## 2 545801 66.41 19.28 14.32 0.05
## 3 545803 66.82 21.68 11.50 0.14
## 4 545805 46.74 41.76 11.50 0.13
## 5 545806 46.74 41.76 11.50 0.13
## 6 545807 46.74 41.76 11.50 0.13# re-create raster attribute table with aggregate soil properties
rat <- merge(ll, p, by.x = 'ID', by.y = 'mukey', sort = FALSE, all.x = TRUE)
# re-pack RAT
levels(mu) <- rat
# convert raster + RAT --> stack of values
s <- deratify(mu, att = c("sandtotal_r", "silttotal_r", "claytotal_r", "awc_r"))
# graphical check
levelplot(
s[[1:3]],
main = 'Sand, Silt, Clay (RV) 0-200cm\nWeighted Average',
margin = FALSE,
scales = list(draw = FALSE),
col.regions = viridis
)
levelplot(
s[[4]],
main = 'AWC (RV) 0-200cm\nWeighted Average',
margin = FALSE,
scales = list(draw = FALSE),
col.regions = viridis
)
# convert to a representative soil texture class
txt.lut <- read.csv('http://soilmap2-1.lawr.ucdavis.edu/800m_grids/RAT/texture_2550.csv')
# make a copy
texture_060 <- s[[1]]
# note: soil textures that aren't present are dropped from factor levels
texture_060[] <- ssc_to_texcl(sand = s$sandtotal_r[], clay = s$claytotal_r[])
# extract RAT
rat <- levels(texture_060)[[1]]
# add colors
rat <- merge(rat, txt.lut, by.x = 'VALUE', by.y = 'class', all.x = TRUE, sort = FALSE)
# fix column order
rat <- rat[, c('ID', 'VALUE', 'hex', 'names')]
# re-pack
levels(texture_060) <- rat
# check: ok
cols <- levels(texture_060)[[1]]$hex
levelplot(texture_060, att = 'names', margin = FALSE, col.regions = cols, main = "Soil Textures") 
writeRaster(texture_060, file = 'C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt/ws32/soiltexture060.tif', overwrite = TRUE, options=c("COMPRESS=LZW"))
texture60raster<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt/ws32/soiltexture060.tif")
### this has been modified to use new SSURGO Inputs with deep-shallow depth classes
texture60raster<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt/ws32/ssurgo-soiltype-class.tif")
## reproject to raster mask instead
ws32<- raster("grassexport/cwt_ws32/ws32.tif")
newproj<- "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"
reproj60 <- projectRaster(from = texture60raster, crs = proj4string(ws32), method = 'ngb')
plot(texture60raster, zlim = c(1,5), main = "Texture Raster with modified Depth Classes")
{plot(reproj60, zlim = c(1,5), main = "Reprojected Texture Raster with modified Depth Classes and Watersheds")
plot(st_geometry(x), add = TRUE, col = NA)}
## prepare 10m resolution raster and resample prepared soil data to 10m
resample10<- raster(resolution=c(10,10), crs = proj4string(reproj60), ext=extent(reproj60))
reproj6010 <- resample(reproj60,resample10, method = 'ngb')
{plot(reproj6010, zlim = c(1,5), main = "Reprojected and Resampled Texture Raster with modified Depth Classes and Watersheds")
plot(st_geometry(x), add = TRUE, col = NA)}
## Crop prepared raster to extent of watershed of interest
## croppedtexture60 <- crop(reproj6010,extent(x[which(x$WS=='32'),]))
#for now cropping is completed based on delineated watershed extent from raster instead of watershed extent from watershed polygons
croppedtexture60 <- crop(reproj6010,extent(ws32))
# croppedtexture60 <- projectRaster(from = croppedtexture60, crs = newproj, method = 'ngb')
## Mask prepared raster using watershed polygon
#masked60 <- mask(croppedtexture60,x[which(x$WS=='32'),])
## for now masking is completed based on delineated watershed extent from raster instead of watershed extent from watershed polygons
ws32polygon<- rasterToPolygons(ws32, fun = function(x){x==1}, dissolve = TRUE, digits = 12)
extent(croppedtexture60)<- extent(ws32)
{plot(reproj6010, zlim = c(1,5), main = "Reprojected and Resampled Texture Raster with modified Depth Classes and Selected Watershed")
plot(ws32polygon, add = TRUE, col = NA)}
masked60 <- mask(croppedtexture60,ws32polygon)
plot(masked60, zlim = c(1,5), main = "Reprojected and Resampled Texture Raster with modified Depth Classes and Masked Selected Watershed")
plot(x[which(x$WS=='32'),])
plot(st_geometry(x[which(x$WS=='32'),]))
{plot(reproj60, zlim = c(1,5))
plot(st_geometry(x[which(x$WS=='32'),]), add = TRUE, col = NA)}
{plot(as.factor(croppedtexture60), zlim = c(1,5))
plot(st_geometry(x[which(x$WS=='32'),]), add = TRUE, col = NA)}
{plot(as.factor(masked60), zlim = c(1,5))
plot(st_geometry(x[which(x$WS=='32'),]), add = TRUE, col = NA)}
{plot(masked60, zlim = c(1,5), main = "Polygon Delineation")
plot(st_geometry(x[which(x$WS=='32'),]), add = TRUE, col = NA)}
{plot(masked60, zlim = c(1,5), main = "Watershed Delineation")
plot(ws32polygon, add = TRUE, col = NA)}
Code Prepared By Dylan Beaudette and modified
setwd(system.file("extdata/", package = "RHESSysIOinR"))
ws32<- raster("grassexport/cwt_ws32/ws32.tif")
# RSS: cell values are map unit keys
r <- raster('Raster soil survey/Coweeta_Final_raster.tif')
## for now, convert to UTM z17
r <- projectRaster(r, ws32, method = 'ngb')
# convert to raster + RAT
r <- ratify(r)
# RAT: raster attribute table
rat <- read.dbf('Raster soil survey/Coweeta_Final_raster.tif.vat.dbf', as.is = TRUE)
names(rat) <- tolower(names(rat))
# musym is the national mu symbol
rss.mu <- st_read('Raster soil survey/rss_nc/rss_nc.gdb', layer = 'mapunit')## Reading layer `mapunit' from data source
## `C:\Users\Carlos\Documents\R\win-library\4.0\RHESSysIOinR\extdata\Raster soil survey\rss_nc\rss_nc.gdb'
## using driver `OpenFileGDB'rss.co <- st_read('Raster soil survey/rss_nc/rss_nc.gdb', layer = 'component')## Reading layer `component' from data source
## `C:\Users\Carlos\Documents\R\win-library\4.0\RHESSysIOinR\extdata\Raster soil survey\rss_nc\rss_nc.gdb'
## using driver `OpenFileGDB'rss.hz <- st_read('Raster soil survey/rss_nc/rss_nc.gdb', layer = 'chorizon')## Reading layer `chorizon' from data source
## `C:\Users\Carlos\Documents\R\win-library\4.0\RHESSysIOinR\extdata\Raster soil survey\rss_nc\rss_nc.gdb'
## using driver `OpenFileGDB'## check for missing symbols
# none missing from RAT
setdiff(rat$MUSYM, rss.mu$musym)## NULL# map unit table contains some extra
setdiff(rss.mu$musym, rat$MUSYM)## [1] "2y9mf" "2y9mh" "2y9mj" "2y9mk" "2y9ml" "2y9mm" "2y9mn" "2y9mp" "2y9mq"
## [10] "2y9mr" "2y9ms" "2y9mt" "2y9mv" "2y9mw" "2y9mx" "2y9mg" "2y9pb" "2y9pc"
## [19] "2y9pd" "2y9p9" "2y9nz" "2y9nl" "2y9nm" "2y9nw" "2y9nx" "2zyxq" "2zyxw"
## [28] "2zyxy" "2zyy5" "2zyxz" "2zyy0" "2zyxx" "2zyxt" "2zyy6" "2y9np" "2y9nt"
## [37] "2y9nv" "2zyxr" "2zyy1" "2zyy2" "2zyy3" "2y9nh" "2y9nj" "2y9nk" "2zyxv"
## [46] "2zyyb" "2zyyc" "2zyyd" "2zyyf" "2y9nn" "2y9p1" "2y9p2" "2y9nq" "2y9nr"
## [55] "2zyy4"## subset child tables
rss.mu <- rss.mu[rss.mu$musym %in% rat$musym, ]
rss.co <- rss.co[rss.co$mukey %in% rss.mu$mukey, ]
rss.hz <- rss.hz[rss.hz$cokey %in% rss.co$cokey, ]
.dominantCondition <- function(i, v) {
i <- i[i$compkind != 'Miscellaneous area', ]
if(nrow(i) < 1) {
return(NULL)
}
fm <- as.formula(sprintf("comppct_r ~ %s", v))
a <- aggregate(fm, data = i, FUN = sum, na.rm = TRUE)
idx <- order(a[['comppct_r']], decreasing = TRUE)[1]
res <- data.frame(
mukey = i$mukey[1],
v = a[[v]][idx]
)
names(res) <- c('mukey', v)
return(res)
}
.dominantValue <- function(i, v) {
i <- i[i$compkind != 'Miscellaneous area', ]
if(nrow(i) < 1) {
return(NULL)
}
idx <- order(i[['comppct_r']], decreasing = TRUE)[1]
res <- data.frame(
mukey = i$mukey[1],
v = i[[v]][idx]
)
names(res) <- c('mukey', v)
return(res)
}
getDominantCondition <- function(x, v) {
res <- lapply(
split(rss.co, rss.co$mukey),
.dominantCondition, v = v
)
res <- do.call('rbind', res)
return(res)
}
getDominantValue <- function(x, v) {
res <- lapply(
split(rss.co.aws050, rss.co.aws050$mukey),
.dominantValue, v = v
)
res <- do.call('rbind', res)
return(res)
}
co.taxpartsize <- getDominantCondition(rss.co, v = 'taxpartsize')
co.spc <- rss.hz
depths(co.spc) <- cokey ~ hzdept_r + hzdepb_r
hzdesgnname(co.spc) <- 'hzname'
# at dZ = 1cm, use awc_r directly
co.spc$aws <- co.spc$awc_r * 1
a <- slab(co.spc, cokey ~ aws, slab.structure = c(0, 50), slab.fun = sum)## Note: aqp::slice() will be deprecated in aqp version 2.0
## --> Please consider using the more efficient aqp::dice()head(a)## variable cokey value top bottom contributing_fraction
## 1 aws 3244721:2807027 2.40 0 50 1
## 2 aws 3244721:2807439 4.02 0 50 1
## 3 aws 3244721:2807442 2.50 0 50 1
## 4 aws 3244721:2807445 3.73 0 50 1
## 5 aws 3244722:2760595 2.40 0 50 1
## 6 aws 3244722:2760596 4.02 0 50 1rss.co.aws050 <- merge(rss.co, a, by = 'cokey', sort = FALSE)[, c('mukey', 'cokey', 'compkind', 'comppct_r', 'value')]
# no NA allowed in wt. mean / etc.
rss.co.aws050 <- rss.co.aws050[!is.na(rss.co.aws050$value), ]
co.aws050 <- getDominantValue(rss.co.aws050, v = 'value')
names(co.aws050) <- c('mukey', 'aws050')
mu.subset <- rss.mu[, c('mukey', 'musym', 'mukind')]
agg <- merge(mu.subset, co.taxpartsize, by = 'mukey', sort = FALSE)
agg <- merge(agg, co.aws050, by = 'mukey', sort = FALSE)
rat <- merge(rat, agg, by = 'musym', sort = FALSE)
# attempt to split out component / series names
rat$co.names <- stri_split_fixed(str = rat$muname, pattern = ',', n = 2, simplify = TRUE)[, 1]
# this only works because all component names are single-word names
rat$co.names <- stri_split_fixed(str = rat$co.names, pattern = ' ', n = 2, simplify = TRUE)[, 1]
# RAT management
ll.original <- levels(r)[[1]]
# merge and pack updated RAT
ll <- merge(ll.original, rat, by.x = 'ID', by.y = 'value', all.x = TRUE, sort = FALSE)
levels(r) <- ll
# convert select attributes to new raster objects via RAT + codes
r.taxpartsize <- deratify(r, att = 'taxpartsize')
plot(r.taxpartsize)
ws32polygon<- rasterToPolygons(ws32, fun = function(x){x==1}, dissolve = TRUE, digits = 12)
maskeddsm <- mask(r.taxpartsize,ws32polygon)
plot(maskeddsm, main = "Masked DSM")
Reclassify Soil Maps
Default Texture Values in RHESSys 1 - Clay 2 - Silty Clay 3 - Silty Clay Loam 4 - Sandy Clay 5 - Sandy Clay Loam 6 - Clay Loam 7 - Silt 8 - Silty Loam 9 - Loam 10 - Sand 11 - Loamy Sand 12 - Sandy Loam 13 - Rock 14 - Water
Shallow Files Created and added to defs folder 19 - Shallow Sandy Clay Loam 23 - Shallow Loam 26 - Shallow Sandy Loam
Default Texture Values used by SDA 1 - loamy sand 2 - sandy loam 3 - loam 4 - sandy clay loam 5 - clay loam
Default Texture used by DSM and quick and dirty conversions 1 - loamy > loam 2 - fine-loamy > silty clay loam 3 - coarse-loamy > sandy loam
# original before depth class substitution
#reclasssoil<- c(1,11,2,12,3,9,4,5,5,6,6,0,7,0,8,0,9,0,10,0,11,0,12,0,13,0,14,0)
##this has been modified for quick implementation of shallow-deep classes
reclasssoil<- c(1,12,2,26,3,9,4,23,5,5,6,19)
reclasssoilmat<- matrix(reclasssoil, ncol = 2, byrow = TRUE)
reclasssoilmat## [,1] [,2]
## [1,] 1 12
## [2,] 2 26
## [3,] 3 9
## [4,] 4 23
## [5,] 5 5
## [6,] 6 19reclassreproj<- reclassify(masked60, reclasssoilmat)
plot(reclassreproj, main = "SSURGO SDA Top 50cm DC Texture Method", zlim = c(1,14), col = viridis(14))
writeRaster(reclassreproj, file = 'grassexport/cwt_ws32/ssurgosoil.tif', overwrite = TRUE, options=c("COMPRESS=LZW"))
#plot(SoilsMap2, main = "SSURGO USACE Texture Method")
reclassdsm<- c(1,9,2,3,3,12,4,0,5,0,6,0,7,0,8,0,9,0,10,0,11,0,12,0,13,0,14,0)
reclassdsmmat<- matrix(reclassdsm, ncol = 2, byrow = TRUE)
reclassdsmmat## [,1] [,2]
## [1,] 1 9
## [2,] 2 3
## [3,] 3 12
## [4,] 4 0
## [5,] 5 0
## [6,] 6 0
## [7,] 7 0
## [8,] 8 0
## [9,] 9 0
## [10,] 10 0
## [11,] 11 0
## [12,] 12 0
## [13,] 13 0
## [14,] 14 0reclassmaskeddsm <- reclassify(maskeddsm, reclassdsmmat)
plot(reclassmaskeddsm, main = "DSM Top 50cm DC Texture Method", zlim = c(1,14), col = viridis(14))
#plot(maskeddsm)
reclassstatic<- c(1,9,2,9,3,9,4,9,5,9,6,9,7,9,8,9,9,9,10,9,11,9,12,9,13,9,14,9)
reclassstaticmat<- matrix(reclassstatic, ncol = 2, byrow = TRUE)
reclassstaticmat## [,1] [,2]
## [1,] 1 9
## [2,] 2 9
## [3,] 3 9
## [4,] 4 9
## [5,] 5 9
## [6,] 6 9
## [7,] 7 9
## [8,] 8 9
## [9,] 9 9
## [10,] 10 9
## [11,] 11 9
## [12,] 12 9
## [13,] 13 9
## [14,] 14 9reclassmaskedstatic <- reclassify(maskeddsm, reclassstaticmat)
plot(reclassmaskedstatic, main = "Static Loam Texture", zlim = c(1,14), col = viridis(14))
NLCD Land Cover Classification Legend
11 Open Water 12 Perennial Ice/ Snow 21 Developed, Open Space 22 Developed, Low Intensity 23 Developed, Medium Intensity 24 Developed, High Intensity 31 Barren Land (Rock/Sand/Clay) 41 Deciduous Forest 42 Evergreen Forest 43 Mixed Forest 51 Dwarf Scrub 52 Shrub/Scrub 71 Grassland/Herbaceous 72 Sedge/Herbaceous 73 Lichens 74 Moss 81 Pasture/Hay 82 Cultivated Crops 90 Woody Wetlands 95 Emergent Herbaceous Wetlands
- Alaska only
Because default RHESSys definition files are being used we will reclassify some of this land cover classification for our landscape
41 Deciduous Forest <- 3 veg_deciduous 42 Evergreen Forest <- 4 veg_evergreen
LANDCOVER<- get_nlcd(template = x.buff.gcs, label = 'coweeta', dataset = "landcover")
#CANOPY<- get_nlcd(template = x.buff.gcs, label = 'coweeta', year = 2016, dataset = "canopy")
IMPERVIOUS<- get_nlcd(template = x.buff.gcs, label = 'coweeta', dataset = "impervious")
plot(LANDCOVER, main = "NLCD Landcover")
#plot(CANOPY)
plot(IMPERVIOUS, main = "NLCD Impervious")
LANDCOVERREPROJ <- projectRaster(from = LANDCOVER, crs = proj4string(ws32polygon), method = 'ngb')
#CANOPYREPROJ <- projectRaster(from = CANOPY, crs = proj4string(ws32polygon), method = 'ngb')
IMPERVIOUSREPROJ <- projectRaster(from = IMPERVIOUS, crs = proj4string(ws32polygon), method = 'ngb')
landcovercrop <- crop(LANDCOVERREPROJ,ws32polygon)
#canopycrop <- crop(CANOPYREPROJ,ws32polygon)
imperviouscrop <- crop(IMPERVIOUSREPROJ,ws32polygon)
maskedlandcover <- mask(landcovercrop,ws32polygon)
#maskedcanopy <- mask(canopycrop,ws32polygon)
maskedimpervious <- mask(imperviouscrop,ws32polygon)
plot(maskedlandcover, main = "Masked Landcover")
#plot(maskedcanopy)
plot(maskedimpervious, main = "Masked Impervious")
{plot(maskedlandcover, main = "Masked Landcover with WS Polygon")
plot(ws32polygon, add = TRUE, col = NA)}
{plot(maskedimpervious, main = "Masked Impervious with WS Polygon")
plot(ws32polygon, add = TRUE, col = NA, border = "white")}
extent(LANDCOVER)## class : Extent
## xmin : 1128045
## xmax : 1136025
## ymin : 1402515
## ymax : 1411275extent(LANDCOVERREPROJ)## class : Extent
## xmin : 271758.1
## xmax : 281385.9
## ymin : 3875918
## ymax : 3886032extent(maskedlandcover)## class : Extent
## xmin : 275077
## xmax : 276123.5
## ymin : 3880946
## ymax : 3881563extent(ws32polygon)## class : Extent
## xmin : 275066.7
## xmax : 276126.7
## ymin : 3880940
## ymax : 3881570# template_read(template)
# would be useful to modify or read RHESSys template in R, currently template_read does not display nicely
# Create Templates
{templatetext<- "1
../defs/basin.def
1
../defs/hillslope.def
1
../defs/zone.def
14
../defs/soil_clay.def
../defs/soil_clayloam.def
../defs/soil_loam.def
../defs/soil_loamysand.def
../defs/soil_rock.def
../defs/soil_sand.def
../defs/soil_sandyclay.def
../defs/soil_sandyclayloam.def
../defs/soil_sandyloam.def
../defs/soil_silt.def
../defs/soil_siltyclay.def
../defs/soil_siltyclayloam.def
../defs/soil_siltyloam.def
../defs/soil_water.def
../defs/soil_shallowloam.def
../defs/soil_shallowsandyclayloam
../defs/soil_shallowsandyloam
1
../defs/lu_undev.def
1
../defs/veg_evergreen.def
1
../clim/cwtws32daymet.base
_world basin_ws32.tif 1
_basin basin_ws32.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
basin_parm_ID dvalue 1
latitude value 35.04
basin_n_basestations dvalue 0
_hillslope subbasins.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
hill_parm_ID dvalue 1
gw.storage value 0.0
gw.NO3 value 0.0
hillslope_n_basestations dvalue 0
_zone patch.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
zone_parm_ID dvalue 1
area area
slope aver slope10.tif
aspect spavg aspect10.tif slope10.tif
precip_lapse_rate value 1.0
e_horizon eqn 0.001 0 e_horizon_1000.tif
w_horizon eqn 0.001 0 w_horizon_1000.tif
zone_n_basestations dvalue 1
zone_basestation_ID dvalue 101
_patch patch.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
soil_parm_ID mode ssurgosoil.tif
landuse_parm_ID dvalue 1
fire_parm_ID value 1.0
area area
slope aver slope10.tif
lna eqn 0.001 0 topidx10_100.tif
Ksat_vertical value 1.0
mpar value 0.12
rz_storage value 0.0
unsat_storage value 0.0
sat_deficit value 0.0
snowpack.water_equivalent_depth value 0.28
snowpack.water_depth value 0.0
snowpack.T value -10.0
snowpack.surface_age value 0.0
snowpack.energy_deficit value -0.5
litter.cover_fraction value 1.0
litter.rain_stored value 0.00001544
litter_cs.litr1c value 0.00001945
litter_ns.litr1n value 0.00000082
litter_cs.litr2c value 0.00316727
litter_cs.litr3c value 0.00021824
litter_cs.litr4c value 0.00532178
soil_cs.soil1c value 0.00238081
soil_ns.sminn value 0.00006014
soil_ns.nitrate value 0.00022041
soil_cs.soil2c value 0.02335700
soil_cs.soil3c value 0.06099644
soil_cs.soil4c value 0.37339245
patch_n_basestations dvalue 0
_canopy_strata patch.tif 1
veg_parm_ID dvalue 3
cover_fraction value 1.0
gap_fraction value 0.0
rootzone.depth value 1.0
snow_stored value 0.0
rain_stored value 0.0
cs.pool value 0.0
cs.leafc value 0.0
cs.dead_leafc value 0.0
cs.leafc_store value 0.0
cs.leafc_transfer value 0.0
cs.live_stemc value 0.0
cs.livestemc_store value 0.0
cs.live.stemc_transfer value 0.0
cs.dead_stem value 0.0
cs.deadstemc_store value 0.0
cs.deadstemc_transfer value 0.0
cs.live_crootc value 0.0
cs.livecrootc_store value 0.0
cs.livecrootc_transfer value 0.0
cs.dead_crootc value 2.0
cs.deadcrootc_store value 0.0
cs.deadcrootc_transfer value 0.0
cs.frootc value 0.0
cs.frootc_store value 0.0
cs.frootc_transfer value 0.0
cs_cwdc value 0.0
epv.prev_leafcalloc value 0.0
ns.pool value 0.0
ns.leafn value 0.0
ns.dead_leafn value 0.0
ns.leafn_store value 0.0
ns.leafn_transfer value 0.0
ns.live_stemn value 0.0
ns.livestemn_store value 0.0
ns.livestemn_transfer value 0.0
ns.dead_stem value 0.0
ns.deadstemn_store value 0.0
ns.deadstemn_transfer value 0.0
ns.live_crootn value 0.0
ns.livecrootn_store value 0.0
ns.livecrootn_transfer value 0.0
ns.dead_crootn value 0.0
ns.deadcrootn_store value 0.0
ns.deadcrootn_transfer value 0.0
ns.frootn value 0.0
ns.frootn_store value 0.0
ns.frootn_transfer value 0.0
ns.cwdn value 0.0
ns.retransn value 0.0
epv.wstress_days dvalue 0
epv.min_fparabs value 0.0
epv.min_vwc value 0.0
canopy_strata_n_basestations dvalue 0
"}
write(templatetext,file = "templates/RHESSYSinR_Coweeta_WS32_Spinup.template")Locate RHESSys Spatial Inputs and Specify Template to be used
Create RHESSys Worldfile and Flowtable for 7.4 based on Spatial Inputs and Template
Raster files and template information are used to setup RHESSys worldfile. Flowtable created by “RHESSysPreprocess” must be converted to work. Unsure why. Preprocessing is controlled by template file.
Template file ‘CWWS32_STATIC.template’ is used to run model normally template file ‘CWWS32_STATICspinup.template’ is used to spin up model. Only difference is artificially extended climate dataset.
#setwd(system.file("extdata/", package = "RHESSysIOinR"))
type = "raster"
typepars = "grassexport/cwt_ws32"
template = "templates/RHESSYSinR_Coweeta_WS32_Spinup.template"
# Suffixes of .world, .flow, and .meta will be appended to the worldfile, flowtable, and metadata files respectively.
name = "CWWS32"
overwrite = TRUE
streams = "streams.tif"
# Optional Flowtable Spatial Data
roads = "roads.tif"
# impervious = "impervious_map"
# roofs = "roofs_map"
header = TRUE
# With certain configurations roads does not work. Returns message: road width cannot be 0 during preprocessing
RHESSysPreprocess(
template = template,
name = name,
type = type,
typepars = typepars,
streams = streams,
overwrite = overwrite,
header = header,
parallel = FALSE,
make_stream = TRUE)## [1] Begin world_gen.R
## Reading in maps
## Finished reading in maps
## [1] Writing worldfile
##
|
| | 0%
|
|======================= | 33%
|
|=============================================== | 67%
|
|======================================================================| 100%
## [1] Created worldfile: ./CWWS32.world
## [1] Created header file: ./CWWS32.hdr
## [1] Begin CreateFlownet.R
## Reading in maps
## Finished reading in maps
## [1] Building flowtable
## [1] Generating unique IDs
## [1] Finding patch neighbors
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] Buildling flowtable list
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] Filling 1 pits.
##
|
| | 0%
|
|======================================================================| 100%
## [1] Writing flowtable
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] Created flowtable: ./CWWS32.flow## Make sure to have parallel turned off when preprocessing, convert flowtables after creation for compatibility
convert_flowtable("CWWS32.flow","CWWS321.flow",make_stream=4)## [1] "Reading in flow table"
##
|
| | 0%
|
|======================= | 33%
|
|=============================================== | 67%
|
|======================================================================| 100%Setup RHESSys Inputs/Outputs for Spinup Note that flowtable is converted flowtable. Using original flow table will result in errors.
Original spinup 0219 to 2017. Modified output to 1776-2017 to test code without very long spin-up. Climate files have been modified to original DAYMET timeseries.
## Read climate file to see length available climate data
read_clim("clim/cwtws32extended.rain", dates_out=TRUE)## [1] "1769 12 27 01" "2021 12 31 24"read_clim("clim/cwtws32daymet.rain", dates_out=TRUE)## [1] "1980 01 01 01" "2021 12 31 24"# Begin RHESSys Spinup
##Short Spinup, original spinup run 1776-2018
input_rhessys = IOin_rhessys_input(
version = rh_path,
tec_file = "tecfiles/tec_daily",
world_file = "CWWS32.world",
world_hdr_prefix = "CWWS32",
flowtable = "CWWS321.flow",
start = "1980 11 1 1",
end = "2018 11 1 1",
output_folder = "out",
output_prefix = "cwws32",
commandline_options = c("-b -g"))
input_tec_data = IOin_tec_std(start = "1981 11 1 1",
end = "2018 11 1 1",
output_state = TRUE)
input_hdr = IOin_hdr(
basin = "defs/basin.def",
hillslope = "defs/hillslope.def",
zone = "defs/zone.def",
soil = c("defs/soil_clay.def","defs/soil_clayloam.def","defs/soil_loam.def","defs/soil_loamysand.def","defs/soil_rock.def","defs/soil_sand.def","defs/soil_sandyclay.def","defs/soil_sandyclayloam.def","defs/soil_sandyloam.def","defs/soil_silt.def","defs/soil_siltyclay.def","defs/soil_siltyclayloam.def","defs/soil_siltyloam.def","defs/soil_water.def", "defs/soil_shallowloam.def", "defs/soil_shallowsandyclayloam.def", "defs/soil_shallowsandyloam.def"),
landuse = "defs/lu_undev.def",
stratum = "defs/veg_deciduous.def",
basestations = "clim/cwtws32daymet.base"
)
##suggested standard hydrologic spinup parameters from rhessys wiki, gw values modified based on previous results from Coweeta Landscape
stdpars<- IOin_std_pars(
m = 1,
k = 10,
soil_dep= 0.4,
m_v = 1,
k_v = 1,
gw1 = 0.15,
gw2 = 0.35)Run Spinup
Spinup runs model on one processor for x number of years. Takes a very long time.
Spinup has been turned off for the purposes of this demonstration. Using previously spun up file ~ 12 hours, run multiple times.
Read Spinup
#spinupoutput<- readin_rhessys_output("out/ws32/cwws32_runspinup")
spinupoutput<- readin_rhessys_output("out/cwws32_runspinup")
spinupoutput$bd$date <- make_date(year=spinupoutput$bd$year, month = spinupoutput$bd$month, day = spinupoutput$bd$day)
spinupoutput$bdg$date <- make_date(year=spinupoutput$bdg$year, month = spinupoutput$bdg$month, day = spinupoutput$bdg$day)
range(spinupoutput$bd$year)## [1] 1981 2018plot(spinupoutput$bd$streamflow, type = "l")
plot(spinupoutput$bd$date,spinupoutput$bd$streamflow, type = "l", col = 'black', main = "Predicted Streamflow",
xlab = "Date", ylab = "Streamflow mm")
plot(spinupoutput$bdg$date, spinupoutput$bdg$lai, type = "l", main = "LAI", xlab = "Date")
{plot(spinupoutput$bdg$soilc~spinupoutput$bdg$date,type = "l", main = "Spinup Soil Carbon", xlab = "Date")
soilctrend <- glm(spinupoutput$bdg$soilc~spinupoutput$bdg$date)
abline(soilctrend, col = 'red')}
{plot(spinupoutput$bdg$date,spinupoutput$bdg$soiln, type = "l", main = "Spinup Soil Nitrogen", xlab = "Date")
soilntrend <- glm(spinupoutput$bdg$soiln~spinupoutput$bdg$date)
abline(soilntrend , col = 'green')}
#prepped<- readin_rhessys_output("G:/Coweeta W32/output/basicspinupws32")Configure RHESSys Inputs/Outputs for a single run.
Use World Statefile created during spin-up. To run model one time at patch scale.
COMMAND LINE OPTIONS
-b Basin output option. Print out response variables for specified basins. -c Canopy stratum output option. Print out response variables for specified strata. -g Grow option. Try to read in dynamic bgc input data and output dynamic bgc parameters. -h Hillslope output option. Print out response variables for specified hillslopes. -p Patch output option. Print out response variables for specified patches. -r Routing option. Gives name of flow_table to define explicit routing connectivity. Also trigger use of explicit routing over TOPMODEL approach. -c Stratum output option. Print out response variables for specified strata.
read_clim("clim/cwtws32daymet.rain", dates_out = TRUE)## [1] "1980 01 01 01" "2021 12 31 24"input_rhessys = IOin_rhessys_input(
version = rh_path,
tec_file = "tecfiles/tec_daily",
world_file = "CWWS32.world.Y2018M10D31H1.state",
world_hdr_prefix = "CWWS32",
flowtable = "CWWS321.flow",
start = "2014 11 1 1",
end = "2018 11 1 1",
output_folder = "out",
output_prefix = "cwws32",
commandline_options = c("-b -g -p"))
# commandline_options = c("-b -g -p 1 128 154914 154914")
## TEC file dictates model output, begin output a year in to allow model SM to stabilize
# do not output_state or worldfile may be overwritten as output is created
input_tec_data = IOin_tec_std(start = "2015 11 1 1",
end = "2018 11 1 1",
output_state = FALSE)
input_hdr = IOin_hdr(
basin = "defs/basin.def",
hillslope = "defs/hillslope.def",
zone = "defs/zone.def",
soil = c("defs/soil_clay.def","defs/soil_clayloam.def","defs/soil_loam.def","defs/soil_loamysand.def","defs/soil_rock.def","defs/soil_sand.def","defs/soil_sandyclay.def","defs/soil_sandyclayloam.def","defs/soil_sandyloam.def","defs/soil_silt.def","defs/soil_siltyclay.def","defs/soil_siltyclayloam.def","defs/soil_siltyloam.def","defs/soil_water.def", "defs/soil_shallowloam.def", "defs/soil_shallowsandyclayloam.def", "defs/soil_shallowsandyloam.def"),
landuse = "defs/lu_undev.def",
stratum = "defs/veg_deciduous.def",
basestations = "clim/cwtws32daymet.base")
## Calibrated Values
#stdpars<- IOin_std_pars(m = 1.12, k = 25.23, soil_dep= 0.35, m_v = 4.07, k_v = 0.59, gw1 = 0.09, gw2 = 0.20)
stdpars<- IOin_std_pars(
m = 1.66,
k = 26.16,
soil_dep= 0.33,
m_v = 0.31,
k_v = 11.40,
gw1 = 0.04,
gw2 = 0.21)Run RHESSys once
Read RHESSys Single Run Output
getwd()## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"singleruntest<- readin_rhessys_output("out/cwws32_runsingle")
plot(singleruntest$bd$streamflow~singleruntest$bd$date, type = "l", main = "Streamflow", xlab = "Date", ylab = "Streamflow", col = 'DarkBlue')
plot(as.Date(singleruntest$bd$date),singleruntest$bd$rz_storage, type = "l", col = 'black', main = "Basin Scale Root Zone Storage",
xlab = "Date", ylab = "mm")
#basin scale soil moisture
plot(singleruntest$bd$rz_storage/singleruntest$bdg$root_depth~singleruntest$bd$date, type = "l", main = "RZ_Storage/Root_Depth x Date", xlab = "Date", ylab = "rz_Storage/root_depth", col = 'brown')
plot(singleruntest$bd$lai~singleruntest$bd$date, type = "l", main = "LAI", xlab = "Date", ylab = "LAI", col = 'DarkGreen')
plot(singleruntest$bd$pet~singleruntest$bd$date, type = "l", main = "Potential Evapotranspiration", xlab = "Date", ylab = "PET", col = 'DarkBlue')
plot(singleruntest$bd$et~singleruntest$bd$date, type = "l", main = "Evapotranspiration", xlab = "Date", ylab = "ET", col = 'darkslategray')
plot((singleruntest$bd$unsat_stor/singleruntest$bd$sat_def_z)~singleruntest$bd$date, type = "l", main = "unsat_stor/sat_def_z", xlab = "Date", ylab = "vwc", col = 'DarkGreen')
Plot and Filter Observed Soil Moisture
## Data has not been pre-prepared and must be cleaned
smws32_1<- read.delim("C:/Users/Carlos/Desktop/ORISE/Soil Moisture Data/CW/knb-lter-cwt.1308.19/CWT_132_1_0.txt", header = T)
smws32_2<- read.delim("C:/Users/Carlos/Desktop/ORISE/Soil Moisture Data/CW/knb-lter-cwt.1308.19/CWT_232_1_1308.txt", header = T)
smws32_3<- read.delim("C:/Users/Carlos/Desktop/ORISE/Soil Moisture Data/CW/knb-lter-cwt.1308.19/CWT_332_1_1308.txt", header = T)
smws32_1$Date<- as.Date(smws32_1$Date)
smws32_2$Date<- as.Date(smws32_2$Date)
smws32_3$Date<- as.Date(smws32_3$Date)
smws32_1 <- smws32_1[,!sapply(smws32_1,is.character)]
smws32_2 <- smws32_2[,!sapply(smws32_2,is.character)]
smws32_3 <- smws32_3[,!sapply(smws32_3,is.character)]
#range(smws32_3$smois30A)
smws32_1mean<- aggregate(smws32_1, by = list(smws32_1$Date), FUN = function(x) mean(x, na.rm = TRUE))
smws32_2mean<- aggregate(smws32_2, by = list(smws32_2$Date), FUN = function(x) mean(x, na.rm = TRUE))
smws32_3mean<- aggregate(smws32_3, by = list(smws32_3$Date), FUN = function(x) mean(x, na.rm = TRUE))
smws32_1mean <- smws32_1mean[order(smws32_1mean$Date),]
smws32_2mean <- smws32_2mean[order(smws32_2mean$Date),]
smws32_3mean <- smws32_3mean[order(smws32_3mean$Date),]
#Filter Soil Moisture Datasets
smws32_1mean <- smws32_1mean[smws32_1mean$smois30A>'0',]
smws32_3mean <- smws32_3mean[smws32_3mean$smois30A<='0.45',]
smws32_3mean <- smws32_3mean[smws32_3mean$Date<='2018-05-15',]
plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", main = "Soil Moisture 30cm Site 1", xlab = "Date", ylab = "Soil Moisture 30A")
plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "blue", main = "Soil Moisture 30cm Site 2", xlab = "Date", ylab = "Soil Moisture 30A")
plot(smws32_3mean$Group.1,smws32_3mean$smois30A, type = "l", col = "red", main = "Soil Moisture 30cm Site 3", xlab = "Date", ylab = "Soil Moisture 30A")
{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Observed Soil Moisture WS32 30cm")
lines(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_3mean$Group.1,smws32_3mean$smois30A, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM1","SM2","SM3"), col=c("red","blue","black"), lty=1)}
smws32_1mean <- smws32_1mean[smws32_1mean$Date>='2017-01-01',]
smws32_2mean <- smws32_2mean[smws32_2mean$Date>='2017-01-01',]
smws32_3mean <- smws32_3mean[smws32_3mean$Date>='2017-01-01',]
smws32_1mean <- smws32_1mean[smws32_1mean$Date<='2018-05-15',]
smws32_2mean <- smws32_2mean[smws32_2mean$Date<='2018-05-15',]
smws32_3mean <- smws32_3mean[smws32_3mean$Date<='2018-05-15',]
{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 30cm")
lines(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_3mean$Group.1,smws32_3mean$smois30A, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM1","SM2","SM3"), col=c("red","blue","black"), lty=1)}
{plot(density(smws32_1mean$smois30A), ylim = c(0,14), xlim = c(0,0.35), col = 'red', main = "Density of Observed Soil Moisture WS32 30cm")
lines(density(smws32_2mean$smois30A), col = 'blue')
lines(density(smws32_3mean$smois30A), col = 'black')
legend("topright",inset = 0.05, legend=c("SM1","SM2","SM3"), col=c("red","blue","black"), lty=1)}
{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 Station 1")
lines(smws32_1mean$Group.1,smws32_1mean$smois30B, type = "l", col = "orange", ylim = c(0,0.4))
lines(smws32_1mean$Group.1,smws32_1mean$smois60A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_1mean$Group.1,smws32_1mean$smois60B, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM1-30A","SM1-30B","SM1-60A","SM1-60B"), col=c("red","orange","blue","black"), lty=1)}
{plot(density(smws32_1mean$smois30A), ylim = c(0,14), xlim = c(0,0.35), col = 'red', main = "Density of Observed Soil Moisture WS32 Station 1")
lines(density(smws32_1mean$smois30B), col = 'orange')
lines(density(smws32_1mean$smois60A), col = 'blue')
lines(density(smws32_1mean$smois60B), col = 'black')
legend("topright",inset = 0.05, legend=c("SM1-30A","SM1-30B","SM1-60A","SM1-60B"), col=c("red","orange","blue","black"), lty=1)}
{plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 Station 2")
lines(smws32_2mean$Group.1,smws32_2mean$smois30B, type = "l", col = "orange", ylim = c(0,0.4))
lines(smws32_2mean$Group.1,smws32_2mean$smois60A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_2mean$Group.1,smws32_2mean$smois60B, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM2-30A","SM2-30B","SM2-60A","SM2-60B"), col=c("red","orange","blue","black"), lty=1)}
##extra filtering for smws32_60a
smws32_2mean <- smws32_2mean %>%
mutate(smois60A = replace(smois60A,Date >= "2017-11-26", NA ))
{plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 Station 2 Filtered")
lines(smws32_2mean$Group.1,smws32_2mean$smois30B, type = "l", col = "orange", ylim = c(0,0.4))
lines(smws32_2mean$Group.1,smws32_2mean$smois60A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_2mean$Group.1,smws32_2mean$smois60B, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM2-30A","SM2-30B","SM2-60A","SM2-60B"), col=c("red","orange","blue","black"), lty=1)}
{plot(smws32_3mean$Group.1,smws32_3mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 Station 3")
lines(smws32_3mean$Group.1,smws32_3mean$smois30B, type = "l", col = "orange", ylim = c(0,0.4))
lines(smws32_3mean$Group.1,smws32_3mean$smois60A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_3mean$Group.1,smws32_3mean$smois60B, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM3-30A","SM3-30B","SM3-60A","SM3-60B"), col=c("red","orange","blue","black"), lty=1)}
soilmoisture30a<- data.frame(smws32_1mean$Date,smws32_1mean$smois30A, smws32_1mean$smois30B,smws32_1mean$smois60A,smws32_1mean$smois60B)
soilmoisture30b<- data.frame(smws32_2mean$Date,smws32_2mean$smois30A, smws32_2mean$smois30B,smws32_2mean$smois60A,smws32_2mean$smois60B)
soilmoisture30c<- data.frame(smws32_3mean$Date,smws32_3mean$smois30A,smws32_3mean$smois30B,smws32_3mean$smois60A,smws32_3mean$smois60B)
merged30a<- merge(soilmoisture30a,soilmoisture30b, by.x = "smws32_1mean.Date", by.y = "smws32_2mean.Date", all = FALSE)
merged30b<- merge(merged30a,soilmoisture30c, by.x = "smws32_1mean.Date", by.y = "smws32_3mean.Date", all = TRUE)
merged30c<- as.data.frame(rowMeans((merged30b[-1]), na.rm= TRUE))
mergedsm<- cbind(merged30b,merged30c)
colnames(mergedsm)<- c("Date","smois30site1a","smois30site1b","smois60site1a","smois60site1b","smois30site2a","smois30site2b","smois60site2a","smois60site2b","smois30site3a","smois30site3b","smois60site3a","smois60site3b","mergedsoilmoisture")
plot(mergedsm$Date, mergedsm$mergedsoilmoisture, type = "l", main = "Average WS32 Soil Moisture", xlab = "Date", ylab = "Merged Soil Moisture All Data")
plot(mergedsm)
## Calculate Calibration/Validation Split based on SM data
# range(mergedsm$Date)
daysofsmdata<- as.numeric(range(mergedsm$Date)[2] - range(mergedsm$Date)[1])
#70/30 split
Caldays<- round(0.7*daysofsmdata)
Valdays<- daysofsmdata-Caldays
Caldates<- as.Date(1)
Caldates[1]<- range(mergedsm$Date)[1]
Caldates[2]<- range(mergedsm$Date)[1]+Caldays
Valdates<- as.Date(1)
Valdates[1]<- range(Caldates)[2]+1
Valdates[2]<- range(mergedsm$Date)[2]Plot and Filter Observed Streamflow Data
Original Datasets are recorded as Water Year, We convert to calendar year
## import daily flow dataset
obsws32<- read.csv('C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32.csv', header = F)
# dataset is recorded as water year, here we convert from water year and note that the water year for this landscape begins in November
## calendar year seems to be correct
#colnames(obsws32)<- c("WYR","Month","Day","flow")
obsws32$WYR <- as.numeric(obsws32$V1)
obsws32$Month <- as.numeric(obsws32$V2)
obsws32$Day <- as.numeric(obsws32$V3)
obsws32$flow <- as.numeric(obsws32$V4)
obsws32$CalendarYear<- ifelse((obsws32$Month=="12"|obsws32$Month=="11"),obsws32$WYR-1,obsws32$WYR)
obsws32$CalendarDate <- as.Date(with(obsws32,paste(CalendarYear,Month,Day,sep="-")),"%Y-%m-%d")
##this removal of 0s may be reason for lag, 0s may be purposefully included in dataset
# instead create timeseries of date from start to end of timeseries and then merge
#obsws32<- obsws32[!(obsws32$flow=="0"),]
obsws32 <- drop_na(obsws32)
plot(obsws32$CalendarDate, obsws32$flow, type = "l", ylab = "Observed Discharge mm/d", xlab = "Date", main = "Observed Discharge WS32 Older Dataset")
obsws32_1<- read.csv('C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2017.csv', header = T)
obsws32_2<- read.csv('C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2018.csv', header = T)
obsws32_3<- read.csv('C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2019.csv', header = T)
obsws32_4<- read.csv('C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2020.csv', header = T)
obsws32_5<- read.csv('C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2021.csv', header = T)
obsws32b<- rbind(obsws32_1,obsws32_2,obsws32_3,obsws32_4,obsws32_5)
obsws32b$CalendarDate <- as.Date(with(obsws32b,paste(year,month,day,sep="-")),"%Y-%m-%d")
plot(obsws32b$CalendarDate,obsws32b$flow, type = "l", xlab = "Date", ylab = "Streamflow", main = "Observed Discharge WS32 Newer Datasets Merged")
dataset1 <- data.frame(Date=obsws32$CalendarDate,flow=obsws32$flow)
dataset2 <- data.frame(Date=obsws32b$CalendarDate,flow=obsws32b$flow)
fullset<- rbind(dataset1,dataset2)
### Streamflow at Coweeta logged as CSM
### conversion from EcohydRology project
## because we are converting from CSM we set watershed area to 1
fullset$discharge_mm<- ConvertFlowUnits(cfs = fullset$flow, WA=1, AREAunits = "mi2")
#{plot(fullset$Date,fullset$flow, type = "l", xlab = "Date", ylab = "Flow", main = "Streamflow WS32")
#lines(fullset$Date,fullset$discharge_mm, col = 'red')}
plot(fullset$Date,fullset$discharge_mm, type = "l", xlab = "Date", ylab = "Discharge mm/d", main = "Observed Discharge WS32 All Datasets Merged")
plot(fullset$Date,fullset$discharge_mm, type = "l", xlab = "Date", ylab = " Log Discharge mm/d", main = "Log of Observed Discharge WS32 All Datasets Merged", log = "y")
fullset$observedstreamflow<- fullset$discharge_mm
write.csv(fullset,"C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/ObservedWS32.csv")
obsws32<- fullset
obsws32sm<- merge(obsws32,mergedsm, by = "Date")
#Caldates
#Valdates
obsws32cal <- obsws32[obsws32$Date >= Caldates[1] & obsws32$Date <= Caldates[2], ]
obsws32val <- obsws32[obsws32$Date >= Valdates[1] & obsws32$Date <= Valdates[2], ]
obsws32smcal <- obsws32sm[obsws32sm$Date >= Caldates[1] & obsws32sm$Date <= Caldates[2], ]
obsws32smval <- obsws32sm[obsws32sm$Date >= Valdates[1] & obsws32sm$Date <= Valdates[2], ]Merge Observed and Simulated Data
## merge values instead of subsetting to match up with simulated dates
singleruntest$bd<- merge(singleruntest$bd,obsws32, by.x = "date", by.y = "Date", all = FALSE)
## merge soil moisture values but include all data
singleruntest$data <-merge(singleruntest$bd,obsws32sm, by.x = "date", by.y = "Date", all = TRUE)Plot Observed and Simulated Data, Plot RHESSys Outputs
## Plot Hydrograph comparing modeled and observed streamflow
{hydrograph(input=singleruntest$bd, streamflow=singleruntest$bd$observedstreamflow, streamflow2 = singleruntest$bd$streamflow, timeSeries = singleruntest$bd$date, precip = singleruntest$bd$streamflow,
P.units = "mm", S.units = "mm normalized by basin area", S1.col = 'Blue', S2.col = 'Red')
legend("topleft",inset = .2, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("blue","red"), lty=1:2)}
{plot(singleruntest$data$date,singleruntest$data$observedstreamflow.x, type = "l", xlab = "Date", ylab = "Streamflow", main = "Predicted vs Observed Streamflow")
lines(singleruntest$data$date, singleruntest$data$streamflow, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:1)}
{plot(singleruntest$data$date,singleruntest$data$observedstreamflow.x, type = "l", xlab = "Date", ylab = "Streamflow", main = "Predicted vs Observed Log Streamflow", log = "y")
lines(singleruntest$data$date, singleruntest$data$streamflow, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:1)}
predictedvsobservedstreamlm<- lm(observedstreamflow~streamflow, data = singleruntest$bd)
{plot(singleruntest$bd$observedstreamflow,singleruntest$bd$streamflow, xlab = "Observed", ylab = "Predicted", main = "Predicted vs Observed Streamflow", xlim = c(0,55), ylim = c(0,55))
abline(a=0, b=1, col = 'RED', lty = 2, lwd = 2)
legend("right",inset = 0.01, legend = paste("R2 =",format(summary(predictedvsobservedstreamlm)$r.squared,digits=3)))}
{plot(singleruntest$data$date,(singleruntest$data$rz_storage/singleruntest$data$rootdepth), type = "l", ylim = c(0,0.35), xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Soil Moisture")
lines(singleruntest$data$date, singleruntest$data$mergedsoilmoisture, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Soil Moisture","Predicted Soil Moisture"), col=c("black","red"), lty=1:1)}
predictedvsobservedsoillm<- lm((rz_storage/rootdepth)~mergedsoilmoisture, data = singleruntest$data)
{plot((singleruntest$data$rz_storage/singleruntest$data$rootdepth),singleruntest$data$mergedsoilmoisture, xlab = "Observed", ylab = "Predicted", main = "Predicted vs Observed Basin Scale Soil Moisture", xlim = c(0.05,0.4), ylim = c(0.05,0.4))
abline(a=0, b=1, col = 'RED', lty = 2, lwd = 2)
legend("right",inset = 0.01, legend = paste("R2 =",format(summary(predictedvsobservedsoillm)$r.squared,digits=3)))}
## Plot other model outputs
par(mfrow=c(5,1), mar = c(2,5,4,2))
plot(singleruntest$bd$date,singleruntest$bd$tmin, type = "l", col = 'BLUE', ylim = c(-15,30), ylab = "Temperature", xlab = "Date")
lines(singleruntest$bd$date,singleruntest$bd$tmax, type = "l", col = 'RED')
plot(singleruntest$bd$date,singleruntest$bd$precip, type = "h", ylab = "Precipitation", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$streamflow, type = "l", ylab = "Streamflow", xlab = "Date", col = 'blue')
plot(singleruntest$bd$date,singleruntest$bd$baseflow, type = "l", ylab = "Baseflow", xlab = "Date", col = 'brown')
plot(singleruntest$bd$date,singleruntest$bd$rz_storage, type = "l", ylab = "RZ Storage", xlab = "Date", col = 'dark green')
plot(singleruntest$bd$date,singleruntest$bd$psn ,type = "l", ylab = "PSN", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$et, type = "l", ylab = "Evaoptranspiration", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$pet, type = "l", ylab = "PET", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$gw.Qout, type = "l", ylab = "Q out", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$gw.storage, type = "l", ylab = "GW Storage", xlab = "Date")
par(mfrow=c(1,1))
### will only plot in growth mode , used to check that vegetation is initialized properly ##
plot(singleruntest$bdg$date, singleruntest$bdg$lai, type = "l", main = "LAI", xlab = "Date")
plot(singleruntest$bdg$date, singleruntest$bdg$soilc, type = "l", main = "Soil Carbon", xlab = "Date")
## plot cal/val sets
{plot(singleruntest$data$date,(singleruntest$data$rz_storage/singleruntest$data$rootdepth), type = "l", ylim = c(0,0.35), xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Soil Moisture")
lines(obsws32smcal$Date,obsws32smcal$mergedsoilmoisture, col = "BLUE")
lines(obsws32smval$Date,obsws32smval$mergedsoilmoisture, col = "RED")
legend("topright",inset = 0, legend=c("Predicted Soil Moisture","Observed Calibration Soil Moisture", "Observed Validation Soil Moisture"), col=c("black","blue","red"), lty=c(1,1,1))}
{plot(singleruntest$data$date,singleruntest$data$streamflow, type = "l", xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Streamflow")
lines(obsws32cal$Date,obsws32cal$observedstreamflow, col = "BLUE")
lines(obsws32val$Date,obsws32val$observedstreamflow, col = "RED")
legend("topright",inset = 0, legend=c("Predicted Streamflow","Observed Calibration Streamflow", "Observed Validation Streamflow"), col=c("black","blue","red"), lty=c(1,1,1))}
Testing Precip input vs model precip
Spatial Data Display
setwd(system.file("extdata/", package = "RHESSysIOinR"))
# Original input maps have been cropped before processing in R, without cropping the RHESSYs preparation in the R environment will not function properly, below are the input maps used by R
ws32<- raster("grassexport/cwt_ws32/basin_ws32.tif")
patches<- raster("grassexport/cwt_ws32/patch.tif")
acc10<- raster("grassexport/cwt_ws32/acc10.tif")
aspect10<- raster("grassexport/cwt_ws32/aspect10.tif")
#bt1000<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/b.t1000.tif")
dem10<- raster("grassexport/cwt_ws32/dem10.tif")
dem10f<- raster("grassexport/cwt_ws32/dem10f.tif")
dir10<- raster("grassexport/cwt_ws32/dir10.tif")
drain10<- raster("grassexport/cwt_ws32/drain10.tif")
#ht1000<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/h.t1000.tif")
hillslope<- raster("grassexport/cwt_ws32/hillslope.tif")
roads<- brick("grassexport/cwt_ws32/roads.tif")
#shade10<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/shade10.tif")
slope10<- raster("grassexport/cwt_ws32/slope10.tif")
streams<- raster("grassexport/cwt_ws32/streams.tif")
topidx10<- raster("grassexport/cwt_ws32/topidx10_100.tif")
#xmap<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/xmap.tif")
#ymap<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/ymap.tif")
plot(ws32, main ="ws32")
plot(patches, main = "patches")
plot(acc10, main = "acc10")
plot(aspect10, main = "aspect10")
#levelplot(bt1000, col.regions = brewer.pal(name='Spectral', n = 11),at=seq(0,6500,1))
#plot(bt1000, main = "bt1000")
plot(dem10, main = "dem10")
plot(dem10f, main = "dem10f")
plot(dir10, main = "dir10")
plot(drain10, main = "drain10")
#plot(ht1000, main = "ht1000")
plot(hillslope, main = "hillslope")
plot(roads, main = "roads")
#plot(shade10, main = "shade10")
plot(slope10, main = "slope10")
plot(streams, main = "streams")
plot(topidx10, main = "topidx10")
#plot(xmap, main = "xmap")
#plot(ymap, main = "ymap")ObservedSMLocationsnt <- read_sf("C:/Users/Carlos/Desktop/ORISE/Soil Moisture Data/CW/knb-lter-cwt.1308.19/1308.kml")
ObservedSMLocations <- st_transform(ObservedSMLocationsnt, crs(dem10))
colorramp<- viridis(128)
{plot(dem10, col = colorramp, main = "Coweeta Observed Soil Moisture Locations")
points(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1]), c(ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", lwd = 1)}
ObservedSMLocations$Name[3]## [1] "132 plot center"ObservedSMLocations[[3]][[3]]## POINT (275432.6 3881522)ObservedSMLocations$Name[4]## [1] "232 plot center"ObservedSMLocations[[3]][[4]]## POINT (275578.2 3881429)ObservedSMLocations$Name[5]## [1] "332 plot center"ObservedSMLocations[[3]][[5]]## POINT (275760.9 3881284)matrix(c(ObservedSMLocations$Name[3],ObservedSMLocations$Name[4],ObservedSMLocations$Name[5],ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), nrow = 3, ncol = 3)## [,1] [,2] [,3]
## [1,] "132 plot center" "275432.649744231" "3881522.05364127"
## [2,] "232 plot center" "275578.214010129" "3881428.55011708"
## [3,] "332 plot center" "275760.864807773" "3881284.18476691"## Extract values for single date
displayday <- singleruntest$pd[which(singleruntest$pd$date=="2018-03-21")]
## reclassify all patches to na and then reclassify to model output values
reclass_all<- c(-2147483648, 2147483647,NA)
reclass_allm<- matrix(reclass_all, ncol = 3, byrow = TRUE)
displaydaypatches<- reclassify(patches,reclass_allm)
displayday$calculatedrzsm <- displayday$rz_storage/displayday$root.depth
reclass_displayday <- cbind(displayday$patchID,displayday$calculatedrzsm)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches<- reclassify(patches,reclass_displayday)
## Plot Spatial Data Output for RZ Storage
#plot(displaydaypatches, xlim = c(280900,282500), ylim = c(4870000,4872000))
{plot(mask(displaydaypatches,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 RZ Soil Moisture Prediction for 1 day", line = -2)
points(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1]), c(ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", cex = 2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}
reclass_displayday <- cbind(displayday$patchID,displayday$root.depth)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches2<- reclassify(patches,reclass_displayday)
{plot(mask(displaydaypatches2,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 root.depth Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}
reclass_displayday <- cbind(displayday$patchID,displayday$rz_storage)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches3<- reclassify(patches,reclass_displayday)
{plot(mask(displaydaypatches3,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 rz_storage Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}
reclass_displayday <- cbind(displayday$patchID,displayday$sat_def)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches4<- reclassify(patches,reclass_displayday)
{plot(mask(displaydaypatches4,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 saturation deficit Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}
## Plot Hydrograph
{hydrograph(input=singleruntest$bd, streamflow=singleruntest$bd$observedstreamflow, streamflow2 = singleruntest$bd$streamflow, timeSeries = singleruntest$bd$date, precip = singleruntest$bd$streamflow,
P.units = "mm", S.units = "mm normalized by basin area", S1.col = 'Blue', S2.col = 'Red')
legend("topleft",inset = .2, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("blue","red"), lty=1:2)}
### end show patch map
patchnumbers<- raster::extract(patches,matrix(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), nrow = 3, ncol = 2))
ObserveSoilMoisturePatches <- matrix(c(ObservedSMLocations$Name[3],ObservedSMLocations$Name[4],ObservedSMLocations$Name[5],ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2],patchnumbers), nrow = 3, ncol = 4)
## Extract values for a single patch, select patch ID
ObserveSoilMoisturePatches[1,4]## [1] "136557"displaypatch<- singleruntest$pd[which(singleruntest$pd$patchID==ObserveSoilMoisturePatches[1,4])]
ObserveSoilMoisturePatches[2,4]## [1] "141942"displaypatch2<- singleruntest$pd[which(singleruntest$pd$patchID==ObserveSoilMoisturePatches[2,4])]
ObserveSoilMoisturePatches[3,4]## [1] "149478"displaypatch3<- singleruntest$pd[which(singleruntest$pd$patchID==ObserveSoilMoisturePatches[3,4])]
plot(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", main = 'CW WS27 More Detailed RZ Storage Patch #147322', ylim = c(0,.4))
{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", ylim = c(0,.4), col = 'darkorange', xlim=as.Date(c("2016-10-1","2019-10-1")), main = "CW WS32 More Detailed Soil Moisture vs Patches",
xlab = "Date", ylab = "Soil Moisture mm")
lines(smws32_2mean$Group.1,smws32_2mean$smois30A,type = "l", col = 'darkred')
lines(smws32_3mean$Group.1,smws32_3mean$smois30A,type = "l", col = 'blue')
lines(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'red', lty = 2)
lines(displaypatch2$date,(displaypatch2$rz_storage/displaypatch2$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'green', lty = 2)
lines(displaypatch3$date,(displaypatch3$rz_storage/displaypatch3$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'black', lty = 2)
legend("bottomleft",inset = .01, legend=c("Observed Patch 1","Predicted Patch 1","Observed Patch 2","Predicted Patch 2","Observed Patch 3","Predicted Patch 3"), col=c("darkorange", "red","darkred","green","blue","black"), lty=c(1,2,1,2,1,2))}
{plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", ylim = c(0,0.4), col = 'BLUE', xlim=as.Date(c("2016-10-1","2019-10-1")), main = "CW WS32 More Detailed Soil Moisture vs Patch #154914",
xlab = "Date", ylab = "Soil Moisture mm")
lines(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", main = 'CW WS27 More Detailed RZ Storage Patch #154914', col = 'red')
legend("bottomleft",inset = .01, legend=c("Observed Soil Moisture","Predicted Soil Moisture"), col=c("blue","red"), lty=1:1)}
#
#obsws32cal <- obsws32[obsws32$Date >= Caldates[1] & obsws32$Date <= Caldates[2], ]
#obsws32val <- obsws32[obsws32$Date >= Valdates[1] & obsws32$Date <= Valdates[2], ]
#obsws32smcal <- obsws32sm[obsws32sm$Date >= Caldates[1] & obsws32sm$Date <= Caldates[2], ]
#obsws32smval <- obsws32sm[obsws32sm$Date >= Valdates[1] & obsws32sm$Date <= Valdates[2], ]
head(obsws32cal)## Date flow discharge_mm observedstreamflow
## 22961 2017-01-01 2.30 2.172645 2.172645
## 22962 2017-01-02 3.04 2.871669 2.871669
## 22963 2017-01-03 2.72 2.569388 2.569388
## 22964 2017-01-04 2.11 1.993165 1.993165
## 22965 2017-01-05 1.92 1.813686 1.813686
## 22966 2017-01-06 1.81 1.709777 1.709777head(obsws32val)## Date flow discharge_mm observedstreamflow
## 23311 2017-12-17 2.14 2.021504 2.021504
## 23312 2017-12-18 2.27 2.144306 2.144306
## 23313 2017-12-19 2.23 2.106521 2.106521
## 23314 2017-12-20 3.39 3.202289 3.202289
## 23315 2017-12-21 2.74 2.588281 2.588281
## 23316 2017-12-22 2.60 2.456033 2.456033head(obsws32smcal)## Date flow discharge_mm observedstreamflow smois30site1a smois30site1b
## 1 2017-01-01 2.30 2.172645 2.172645 0.2194375 0.2468125
## 2 2017-01-02 3.04 2.871669 2.871669 0.2241528 0.2509896
## 3 2017-01-03 2.72 2.569388 2.569388 0.2184132 0.2409236
## 4 2017-01-04 2.11 1.993165 1.993165 0.1952743 0.2236250
## 5 2017-01-05 1.92 1.813686 1.813686 0.1825451 0.2149410
## 6 2017-01-06 1.81 1.709777 1.709777 0.1745729 0.2097500
## smois60site1a smois60site1b smois30site2a smois30site2b smois60site2a
## 1 0.1988021 0.2463715 0.2162465 0.2141181 0.2191528
## 2 0.2243646 0.2737917 0.2244410 0.2392847 0.2239028
## 3 0.2177813 0.2657361 0.2209722 0.2291771 0.2154028
## 4 0.1993090 0.2483194 0.2058437 0.2039549 0.1967431
## 5 0.1845556 0.2335799 0.1972882 0.1883021 0.1870486
## 6 0.1744132 0.2231181 0.1914861 0.1774236 0.1808646
## smois60site2b smois30site3a smois30site3b smois60site3a smois60site3b
## 1 0.1775868 0.2830729 0.2573576 0.2625312 0.2399167
## 2 0.2186528 0.2833924 0.2626007 0.2942917 0.2684410
## 3 0.2180486 0.2697604 0.2480833 0.2926944 0.2704132
## 4 0.1968403 0.2399167 0.2206285 0.2794965 0.2609514
## 5 0.1786215 0.2264861 0.2081111 0.2668299 0.2512847
## 6 0.1669340 0.2187014 0.2009896 0.2569028 0.2437292
## mergedsoilmoisture
## 1 0.2317839
## 2 0.2490255
## 3 0.2422839
## 4 0.2225752
## 5 0.2099661
## 6 0.2015738head(obsws32smval)## Date flow discharge_mm observedstreamflow smois30site1a smois30site1b
## 351 2017-12-17 2.14 2.021504 2.021504 0.1181493 0.2395313
## 352 2017-12-18 2.27 2.144306 2.144306 0.1165625 0.2392708
## 353 2017-12-19 2.23 2.106521 2.106521 0.1135139 0.2332326
## 354 2017-12-20 3.39 3.202289 3.202289 0.1484514 0.2609861
## 355 2017-12-21 2.74 2.588281 2.588281 0.1391528 0.2532083
## 356 2017-12-22 2.60 2.456033 2.456033 0.1312049 0.2428194
## smois60site1a smois60site1b smois30site2a smois30site2b smois60site2a
## 351 0.1677882 0.2054826 0.1575382 0.1352778 NA
## 352 0.1635625 0.2064132 0.1605903 0.1370000 NA
## 353 0.1600174 0.2057812 0.1585660 0.1363750 NA
## 354 0.1844479 0.2319132 0.1872778 0.1721458 NA
## 355 0.2092639 0.2517917 0.1925139 0.1736111 NA
## 356 0.1954965 0.2400937 0.1847604 0.1649479 NA
## smois60site2b smois30site3a smois30site3b smois60site3a smois60site3b
## 351 0.09663194 0.2347188 0.1869757 0.2350000 0.2169896
## 352 0.09784375 0.2385694 0.1832847 0.2350000 0.2167847
## 353 0.09828472 0.2289757 0.1792153 0.2350000 0.2160000
## 354 0.12405208 0.2609826 0.2203681 0.2421354 0.2173993
## 355 0.14059028 0.2525278 0.2136979 0.2712326 0.2365937
## 356 0.13011806 0.2397049 0.1995382 0.2683160 0.2404201
## mergedsoilmoisture
## 351 0.1812803
## 352 0.1813529
## 353 0.1786329
## 354 0.2045600
## 355 0.2121985
## 356 0.2034018#subset Validation range based on Val Dates set by available SM data
validationsubset <- singleruntest$bd[singleruntest$bd$date >= Valdates[1] & singleruntest$bd$date <= Valdates[2], ]
## Fit tests for streamflow
NSE(validationsubset$streamflow,validationsubset$observedstreamflow)## [1] 0.5994421NSE(validationsubset$streamflow,validationsubset$observedstreamflow, FUN = log, epsilon = "Pushpalatha2012", na.rm=TRUE)## [1] 0.722282KGE(validationsubset$streamflow,validationsubset$observedstreamflow)## [1] 0.802176plot(validationsubset$streamflow~validationsubset$date, type = "l", col = 'red')
lines(validationsubset$observedstreamflow~validationsubset$date)
# Plot Observed and Predicted Streamflow
{plot(validationsubset$date,validationsubset$streamflow, type = 'l', lty = 2, col = 'RED', main = "Streamflow")
lines(validationsubset$date,validationsubset$observedstreamflow, type = 'l', col = 'BLACK')
legend("topright",inset = .1, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:2)}
{plot(singleruntest$bd$date,singleruntest$bd$streamflow, type = 'l', lty = 2, col = 'RED', main = "Streamflow")
lines(singleruntest$bd$date,singleruntest$bd$observedstreamflow, type = 'l', col = 'BLACK')
legend("topright",inset = .1, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:2)}
#Plot Model RZ SM outputs for patch
{plot(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l",ylim = c(0,0.50), ylab = "Root Zone Storage (mm/m) ", xlab = "Date", main = 'Root Zone Display Patch')
lines(displaypatch$date,(displaypatch$rz_field_capacity/displaypatch$root.depth), type = "l", col = "blue")
lines(displaypatch$date,(displaypatch$rz_wilting_point/displaypatch$root.depth), type = "l", col = "red")
legend("bottomleft",inset = .1, legend=c("Predicted RZ Storage","RZ Field Capacity","RZ Wilting Point"), col=c("black","blue","red"), lty=1:1)}
Setup RHESSys Calibration, CLHS used to reduce amount of samples processed
The model calibrates over the following parameters.
M - Decay of hydraulic conductivity with depth K - Hydraulic conductivity at the surface Sd - soil depth Mv - Vertical decay of hydraulic conductivity with depth Kv - Vertical hydraulic conductivity at the surface GW1 - Saturated to Groundwater Coefficient GW2 - Groundwater Loss Coefficient
Parameter & Typical Range from RHESSys Wiki
M - .01 - 20 K - 1 - 600 GW1 - .001 - .3 GW2 - .01 - .9
" In order to generate an adequate sample across the full range for each parameter, values from .001 - .01, values from .01 - 1, and values over 1 should be generated separately. The m and K parameters are multipliers on the initial values set for m and K (values assigned to the m and K maps)." - Taken from RHESSys Wiki
setwd(system.file("extdata/", package = "RHESSysIOinR"))
# number of samples for full data space
fs <- 10000
# number of samples to calibrate over, 100 samples take approx 40 mins in base R and 1 hour in rmarkdown
n <- 150
# create data frame
d <- data.frame(
m = runif(fs, min=0.5, max=3),
k = runif(fs, min=0.5, max=60),
soil_dep = runif(fs, min=0.1, max=0.5),
m_v = runif(fs, min=0.5, max=15),
k_v = runif(fs, min=0.5, max=15),
gw1 = runif(fs, min=0.1, max=0.3),
gw2 = runif(fs, min=0.1, max=0.5))
# marginal distributions, for analysis
m <- melt(d)## No id variables; using all as measure variablesbwplot(variable ~ value, data = m, par.settings = tactile.theme(), main = "distribution of calibration parameters before clhs subsetting")
# cLHS subset
s <- clhs(d, size = n, simple = TRUE, use.cpp = TRUE)
# combine full dataset + subset
g <- make.groups(full = d, cLHS = d[s, ])
# create new list
params<- d[s, ]
params[1:5,]## m k soil_dep m_v k_v gw1 gw2
## 1152 0.5114497 21.44952 0.1439679 9.6779017 9.5969864 0.1829364 0.4490044
## 1798 0.9612092 16.05240 0.4167406 4.8633875 14.2459499 0.2539121 0.2319639
## 5984 1.3464745 52.34296 0.4217323 1.3377789 14.4170302 0.2133951 0.4848750
## 7019 0.5357427 57.27382 0.1720589 0.5526789 6.8889040 0.1237745 0.3778741
## 8216 2.8337005 51.29474 0.1356276 2.2444764 0.6041038 0.1726981 0.2726171## Quick and dirty way to include "uncalibrated" result to test against, groundwater model still included for now, unsure of what to do with parameters for gw
## params[1,] = c(1,1,1,1,1,0,0)
## did not produce significantly different result
#write param table to specific text file if file does not exist
## if name is not changed then old table is used
# paramtablename = "out/cwparamtable021723.txt"
####WORKING
### WORKING
### WORKING
getwd()## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"write.table(params, file = "out/cwparamtable032723.txt")
### writes new params table with todays date to a text file
#write.table(params, file = paste("out/cwparamtable",format(Sys.Date(),"%m%d%y"),".txt",sep = ""))
## reads in specific params file, used for keeping input params constant for analysis
## do not rewrite file if it exits
#{if (file.exists("out/cwparamtable102522.txt"))
# print("text")}
#params <- read.table(file = "out/cwparamtable102522.txt")
params <- read.table(file = "out/cwparamtable032823.txt")
# prepare RHESSys for run
params<- as.data.frame(params)
input_rhessys = IOin_rhessys_input(
version = rh_path,
tec_file = "tecfiles/tec_daily",
world_file = "CWWS32.world.Y2018M10D31H1.state",
world_hdr_prefix = "CWWS32",
flowtable = "CWWS321.flow",
start = "2014 11 1 1",
end = "2018 11 1 1",
output_folder = "out",
output_prefix = "cwws32",
commandline_options = c("-b -g "))
#commandline_options = c("-b -p 1 128 154914 154914 ")) <- Patch Specific Output
## must incorporate single patch output or local patch output for the purposes of calibration
## make sure to run calibration with no modifiers or 1 for everything except gw1 and gw2
## must refine calibration to do multiple steps, first gw1 gw2, then other variables
## multi point calibration and validation
input_tec_data = IOin_tec_std(start = "2015 11 1 1",
end = "2018 11 1 1",
output_state = FALSE)
input_hdr = IOin_hdr(
basin = "defs/basin.def",
hillslope = "defs/hillslope.def",
zone = "defs/zone.def",
soil = c("defs/soil_clay.def","defs/soil_clayloam.def","defs/soil_loam.def","defs/soil_loamysand.def","defs/soil_rock.def","defs/soil_sand.def","defs/soil_sandyclay.def","defs/soil_sandyclayloam.def","defs/soil_sandyloam.def","defs/soil_silt.def","defs/soil_siltyclay.def","defs/soil_siltyclayloam.def","defs/soil_siltyloam.def","defs/soil_water.def", "defs/soil_shallowloam.def", "defs/soil_shallowsandyclayloam.def", "defs/soil_shallowsandyloam.def"),
landuse = c("defs/lu_undev.def","defs/lu_agriculture.def","defs/lu_raingarden.def","defs/lu_urban.def"),
stratum = c("defs/veg_deciduous.def","defs/veg_deciduous_BES.def","defs/veg_eucalypt.def","defs/veg_evergreen.def","defs/veg_grass.def","defs/veg_lawn_2cm.def","defs/veg_lawn_5cm.def","defs/veg_lawn_10cm.def","defs/veg_nonveg.def"),
basestations = "clim/cwtws32daymet.base")Run RHESSys n times sequentially
Run RHESSys in Parallel using the foreach package, leaving 4 cores available for other tasks
Takes about 1 hour for 100 runs at 10x10m Daily WS32.
Read in RHESSys Calibration Outputs and Plot
Calibration is performed using GLUE methodology and Nash-Sutcliffe Model Efficiency (NSE), Log-Normal Nash-Sutcliffe Model Efficiency (lnNSE), and Kling-Gupta Efficiency (KGE).
setwd(system.file("extdata/", package = "RHESSysIOinR"))
# Read in RHESSys Calibration runs
for(i in 1:n) { assign(paste0("cwws32_run",i), readin_rhessys_output(paste0("out/cwws32_run",i)))}
#Prepare observed Streamflow data
obsws32$Date<- as.Date(obsws32$Date)
## plot observed data for period of record and calibration/validation periods as defined by SM data
{plot(obsws32$Date,obsws32$discharge_mm, type = 'l', xlim=c(as.Date(input_rhessys$start_date, "%Y%m%d"),as.Date(input_rhessys$end_date, "%Y%m%d")),ylab = "Observed Discharge mm/d", xlab = "Date", main = "Observed Streamflow Discharge for period of time model is run")
abline(v=as.Date(Caldates[1]), col = 'red', lwd = 2, lty = 2)
abline(v=as.Date(Caldates[2]), col = 'red', lwd = 2, lty = 2)
abline(v=as.Date(Valdates[1]), col = 'blue', lwd = 2, lty = 5)
abline(v=as.Date(Valdates[2]), col = 'blue', lwd = 2, lty = 5)
legend("topright",inset = 0.05, legend=c("Calibration","Validation"), col=c("red","blue"), lty=1)}
##merge calibration runs with observed so dates match up, instead of subsetting
for(i in 1:n) { assign(paste0("simmerge",i), merge(obsws32,eval(parse(text = paste0("cwws32_run",i,"$bd"))), by.x = "Date", by.y="date", all = FALSE))
assign(paste0("calsubset",i),eval(parse(text = paste0("simmerge",i)))[eval(parse(text = paste0("simmerge",i,"$Date"))) >= Caldates[1] & eval(parse(text = paste0("simmerge",i,"$Date"))) <= Caldates[2], ])}
for(i in 1:n) { assign(paste0("simmergesm",i), merge(obsws32sm,eval(parse(text = paste0("cwws32_run",i,"$bd"))), by.x = "Date", by.y="date", all = FALSE))
assign(paste0("calsubsetsm",i),eval(parse(text = paste0("simmergesm",i)))[eval(parse(text = paste0("simmergesm",i,"$Date"))) >= Caldates[1] & eval(parse(text = paste0("simmergesm",i,"$Date"))) <= Caldates[2], ])}
# simmerge1<- merge(obsws32,eval(parse(text = paste0("cwws32_run",i,"$bd"))), by.x = "Date", by.y="date", all = FALSE)
#simmergesm1<- merge(observedsubsetsm,eval(parse(text = paste0("cwws32_run",i,"$bd"))), by.x = "Date", by.y="date", all = FALSE)
#subset observed data to calibration dates
# calsubset <- eval(parse(text = paste0("simmerge",i)))[eval(parse(text = paste0("simmerge",i,"$Date"))) >= Caldates[1] & eval(parse(text = paste0("simmerge",i,"$Date"))) <= Caldates[2], ]
#calsubsetsm <- observedsubsetsm[observedsubsetsm$Date >= Caldates[1] & observedsubsetsm$Date <= Caldates[2], ]setwd(system.file("extdata/", package = "RHESSysIOinR"))
NSElist<- c()
lnNSElist<- c()
KGElist<- c()
smNSElist<- c()
smlnNSElist<- c()
smKGElist<- c()
for(i in 1:n) { assign(paste0("NSEobs",i), NSE(sim = as.numeric(eval(parse(text = paste0("simmerge",i,"$streamflow")))), obs = as.numeric(eval(parse(text = paste0("simmerge",i,"$observedstreamflow")))) ))
NSElist[[i]]<- eval(parse(text = paste0("NSEobs",i))) }
for(i in 1:n) { assign(paste0("lnNSEobs",i), NSE( sim = as.numeric(eval(parse(text = paste0("simmerge",i,"$streamflow")))), obs = as.numeric(eval(parse(text = paste0("simmerge",i,"$observedstreamflow")))), FUN = log, epsilon = "Pushpalatha2012", na.rm=TRUE))
lnNSElist[[i]]<- eval(parse(text = paste0("lnNSEobs",i))) }
for(i in 1:n) { assign(paste0("KGEobs",i), KGE(sim = as.numeric(eval(parse(text = paste0("simmerge",i,"$streamflow")))), obs = as.numeric(eval(parse(text = paste0("simmerge",i,"$observedstreamflow"))))))
KGElist[[i]]<- eval(parse(text = paste0("KGEobs",i))) }
for(i in 1:n) { assign(paste0("smNSEobs",i), NSE(sim = as.numeric(eval(parse(text = paste0("simmergesm",i,"$rz_storage","/simmergesm",i,"$rootdepth")))), obs = as.numeric(eval(parse(text = paste0("simmergesm",i,"$mergedsoilmoisture"))))))
smNSElist[[i]]<- eval(parse(text = paste0("smNSEobs",i))) }
for(i in 1:n) { assign(paste0("smlnNSEobs",i), NSE(sim = as.numeric(eval(parse(text = paste0("simmergesm",i,"$rz_storage","/simmergesm",i,"$rootdepth")))), obs = as.numeric(eval(parse(text = paste0("simmergesm",i,"$mergedsoilmoisture")))), FUN = log, epsilon = "Pushpalatha2012", na.rm=TRUE))
smlnNSElist[[i]]<- eval(parse(text = paste0("smlnNSEobs",i))) }
for(i in 1:n) { assign(paste0("smKGEobs",i), KGE(sim = as.numeric(eval(parse(text = paste0("simmergesm",i,"$rz_storage","/simmergesm",i,"$rootdepth")))), obs = as.numeric(eval(parse(text = paste0("simmergesm",i,"$mergedsoilmoisture"))))))
smKGElist[[i]]<- eval(parse(text = paste0("smKGEobs",i))) }
NSElist[1:5]## [[1]]
## [1] -0.003875825
##
## [[2]]
## [1] 0.0767704
##
## [[3]]
## [1] 0.1992447
##
## [[4]]
## [1] -0.01049478
##
## [[5]]
## [1] 0.4550887KGElist[1:5]## [[1]]
## [1] 0.4662955
##
## [[2]]
## [1] 0.495591
##
## [[3]]
## [1] 0.5451089
##
## [[4]]
## [1] 0.4611981
##
## [[5]]
## [1] 0.5980345lnNSElist[1:5]## [[1]]
## [1] 0.2647235
##
## [[2]]
## [1] 0.2197867
##
## [[3]]
## [1] 0.2514198
##
## [[4]]
## [1] 0.1662837
##
## [[5]]
## [1] 0.3505744Reduce(max, NSElist)## [1] 0.4876977Reduce(max, lnNSElist)## [1] 0.3570073Reduce(max, KGElist)## [1] 0.6270194Reduce(max, smNSElist)## [1] -1.624272Reduce(max, smlnNSElist)## [1] -1.50897Reduce(max, smKGElist)## [1] 0.002953977which(NSElist == Reduce(max, NSElist), arr.ind=TRUE)## [1] 142which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE)## [1] 58which(KGElist == Reduce(max, KGElist), arr.ind=TRUE)## [1] 65which(smNSElist == Reduce(max, smNSElist), arr.ind=TRUE)## [1] 72which(smlnNSElist == Reduce(max, smlnNSElist), arr.ind=TRUE)## [1] 34which(smKGElist == Reduce(max, smKGElist), arr.ind=TRUE)## [1] 4paramtable<- read.table(file = "out/cwparamtable032823.txt")
paramtable[which(NSElist == Reduce(max, NSElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2
## 4947 1.478595 6.810474 0.3167354 1.225985 14.37205 0.266005 0.1388771paramtable[which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2
## 6586 1.996897 11.9516 0.1659503 14.59814 11.96051 0.04187779 0.3730755paramtable[which(KGElist == Reduce(max, KGElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2
## 5563 1.669432 26.16634 0.3261013 0.3147534 11.40159 0.03675462 0.2058475paramtable[which(smNSElist == Reduce(max, smNSElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2
## 960 0.8611576 1.628167 0.4270699 0.5689802 12.58125 0.2944458 0.2843299paramtable[which(smlnNSElist == Reduce(max, smlnNSElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2
## 6887 1.241347 2.344012 0.4243403 13.31544 6.071543 0.1391299 0.3108362paramtable[which(smKGElist == Reduce(max, smKGElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2
## 9148 0.6367556 6.600494 0.2709346 0.7792454 11.35695 0.262354 0.4609033# paramtable[7,]
###plot all data for GLUE
paramtable$NSE<- NSElist
paramtable$lnNSE<- lnNSElist
paramtable$KGE<- KGElist
paramtable$run<- seq.int(nrow(paramtable))
paramtable$smNSE<- smNSElist
paramtable$smlnNSE<- smlnNSElist
paramtable$smKGE<- smKGElist
paramtable<- as.data.frame(lapply(paramtable,unlist))
top100<- paramtable[order(-paramtable$lnNSE),][1:100,]
top100$run## [1] 58 19 126 125 57 137 40 117 5 76 16 94 42 50 79 128 80 127
## [19] 146 51 49 8 62 142 108 92 135 26 66 41 78 69 87 101 91 111
## [37] 45 73 88 27 71 133 136 144 99 89 21 28 107 67 22 81 110 122
## [55] 20 7 68 85 114 33 145 15 25 6 119 118 53 140 59 98 121 143
## [73] 44 123 120 124 30 116 113 109 139 102 65 47 83 134 14 104 82 29
## [91] 112 100 64 131 77 48 63 84 36 52#use values from top 100 run to table of daily streamflow predictions, for each day predict 95% prediction interval and mean then plot
#subset streamflow value for each day - day 1 day 2 etc
#put all values into a table
top100[["run"]]## [1] 58 19 126 125 57 137 40 117 5 76 16 94 42 50 79 128 80 127
## [19] 146 51 49 8 62 142 108 92 135 26 66 41 78 69 87 101 91 111
## [37] 45 73 88 27 71 133 136 144 99 89 21 28 107 67 22 81 110 122
## [55] 20 7 68 85 114 33 145 15 25 6 119 118 53 140 59 98 121 143
## [73] 44 123 120 124 30 116 113 109 139 102 65 47 83 134 14 104 82 29
## [91] 112 100 64 131 77 48 63 84 36 52toppredictions<- data.frame(cwws32_run1$bd$date)
for (i in top100$run){
# print(i)
toppredictions[,ncol(toppredictions)+1] <- eval(parse(text=paste0("cwws32_run",i,"$bd$streamflow")))
colnames(toppredictions)[ncol(toppredictions)] <- paste0("stream",i)
}
head(toppredictions[1:5,])## cwws32_run1.bd.date stream58 stream19 stream126 stream125 stream57 stream137
## 1 2015-11-01 7.301086 7.621391 7.234245 7.033243 7.219342 6.935841
## 2 2015-11-02 8.339411 8.481753 8.223325 8.076945 8.454141 8.075405
## 3 2015-11-03 7.622130 7.759192 7.172434 7.837778 8.350942 8.026212
## 4 2015-11-04 7.194301 7.304438 6.680168 7.257826 7.691580 7.497977
## 5 2015-11-05 7.160050 7.311759 6.628511 7.158218 7.412052 7.315987
## stream40 stream117 stream5 stream76 stream16 stream94 stream42 stream50
## 1 6.782780 7.173283 6.999768 7.718024 7.229691 7.243884 6.920345 7.543200
## 2 7.845375 8.105602 7.945108 8.758111 8.319140 8.427707 8.165312 8.586470
## 3 7.756190 7.604787 7.850628 7.978340 7.729866 8.074425 8.206078 7.850840
## 4 7.245324 7.235594 7.204286 7.431195 7.370403 7.595946 7.660926 7.434651
## 5 7.036828 7.255528 7.144896 7.330426 7.420635 7.441549 7.419535 7.375747
## stream79 stream128 stream80 stream127 stream146 stream51 stream49 stream8
## 1 6.842465 7.031354 6.967658 6.908412 7.002993 6.982513 7.046788 7.197182
## 2 7.613180 7.979749 7.936525 7.709666 7.807816 7.912854 8.251385 8.231783
## 3 7.243560 7.685994 7.711997 7.469563 7.302992 7.773524 8.112560 7.743575
## 4 6.811250 7.303197 7.208213 6.940977 7.009744 7.280964 7.649492 7.392883
## 5 6.737169 7.224022 7.158457 6.960004 7.127518 7.109254 7.494514 7.438454
## stream62 stream142 stream108 stream92 stream135 stream26 stream66 stream41
## 1 6.890182 7.480707 7.050082 6.994545 8.149884 6.128937 7.076216 7.746484
## 2 7.745624 8.359867 7.697038 8.156216 9.069572 7.083323 8.226110 8.935218
## 3 7.625534 7.664917 7.360530 7.866016 7.740992 6.885493 7.897024 8.100632
## 4 7.266316 7.213749 7.001604 7.509209 7.264247 6.467991 7.464590 7.517883
## 5 7.186921 7.162419 7.042933 7.440129 7.260217 6.335203 7.384197 7.443158
## stream78 stream69 stream87 stream101 stream91 stream111 stream45 stream73
## 1 6.794378 7.119836 7.011195 6.720076 7.547154 6.723873 6.934614 9.046141
## 2 7.660612 7.907789 7.858544 7.854581 8.689763 7.517447 7.604283 10.119782
## 3 7.610500 7.633965 7.695322 7.900355 8.029957 6.980677 7.311125 9.198160
## 4 7.036635 7.211203 7.264788 7.409750 7.537334 6.705402 7.023489 8.598048
## 5 7.061996 7.198216 7.109192 7.301078 7.457512 6.901077 7.142527 8.370346
## stream88 stream27 stream71 stream133 stream136 stream144 stream99 stream89
## 1 6.696809 7.066230 7.179298 6.772535 9.074359 6.780157 6.669808 6.877207
## 2 7.753083 8.544112 8.333969 7.501724 9.963068 7.873735 7.499996 7.675935
## 3 7.602220 8.365502 8.008033 6.912154 7.925216 7.414847 7.176466 7.560846
## 4 7.336603 7.767062 7.515260 6.728857 7.254993 7.135348 6.913612 7.235897
## 5 7.307955 7.521184 7.458419 6.866491 7.204661 7.105699 6.960972 7.206261
## stream21 stream28 stream107 stream67 stream22 stream81 stream110 stream122
## 1 6.905437 6.645647 7.209404 6.768421 6.878164 6.894000 6.484788 7.172463
## 2 8.396006 7.596875 8.336647 7.498912 7.800314 7.756850 7.345957 8.224653
## 3 8.226518 7.715250 7.948836 7.313388 7.270891 7.206889 6.939228 7.731809
## 4 7.674125 7.161561 7.541502 7.062653 6.527299 7.006527 6.764322 7.364885
## 5 7.481645 7.145523 7.480826 7.150338 6.145675 7.169440 6.895039 7.450900
## stream20 stream7 stream68 stream85 stream114 stream33 stream145 stream15
## 1 6.890570 7.051987 7.123813 6.690226 7.013109 7.260308 6.551573 6.539392
## 2 7.904469 8.205911 8.342977 7.342435 7.974642 8.258143 7.489137 7.057671
## 3 7.505268 7.950429 8.033743 7.126135 7.487044 7.591996 7.239523 6.481352
## 4 7.318275 7.511801 7.616806 6.868980 7.244088 7.344860 7.016527 6.276343
## 5 7.425878 7.432365 7.493687 7.006537 7.407971 7.456604 7.122548 6.412394
## stream25 stream6 stream119 stream118 stream53 stream140 stream59 stream98
## 1 6.454897 6.731981 7.288448 6.697253 6.753820 6.710789 7.173417 6.622145
## 2 7.339007 7.870349 8.402888 7.489018 7.329145 7.718049 8.507835 7.299454
## 3 6.869175 8.070802 7.945091 7.234954 7.021659 7.934160 8.002741 6.976581
## 4 6.729247 7.459924 7.570510 6.919776 6.733350 7.424889 7.578688 6.747043
## 5 6.899262 7.264626 7.512899 6.909700 6.882012 7.272107 7.531734 6.941437
## stream121 stream143 stream44 stream123 stream120 stream124 stream30 stream116
## 1 7.066062 7.158570 7.370372 8.067720 6.041898 7.276809 7.248389 6.801399
## 2 8.075929 8.437932 8.183350 9.316698 6.908755 8.509898 8.522186 7.648686
## 3 7.549531 8.219739 7.664192 8.133987 6.222903 7.870076 8.105675 7.320784
## 4 7.329247 7.653680 7.728082 7.512038 6.124654 7.566975 7.698694 6.945106
## 5 7.473999 7.550356 8.447302 7.424370 6.416691 7.614377 7.597161 6.954804
## stream113 stream109 stream139 stream102 stream65 stream47 stream83 stream134
## 1 6.636456 6.606008 6.656626 7.891642 7.290588 7.798762 6.601900 7.746128
## 2 7.359144 7.289312 7.313785 9.290798 7.519508 9.089471 7.074206 8.950239
## 3 7.188012 7.038974 7.106471 8.240449 6.054104 8.014270 6.763529 8.003373
## 4 6.864055 6.837338 6.925474 7.678467 5.469140 7.481455 6.589236 7.595872
## 5 6.867512 6.936015 7.039053 7.584114 5.396264 7.433285 6.791125 7.641679
## stream14 stream104 stream82 stream29 stream112 stream100 stream64 stream131
## 1 8.355276 7.406120 7.747981 7.056639 6.655416 7.285621 6.516596 7.940796
## 2 9.638516 8.769455 8.316809 7.520984 7.387885 8.541844 7.289585 9.125171
## 3 8.153606 8.202823 7.994366 6.672738 7.050887 8.056511 7.409187 7.844993
## 4 7.606671 7.699445 7.583742 6.337906 6.710635 7.659458 6.949297 7.451631
## 5 7.583419 7.620370 7.426955 6.275115 6.741510 7.637652 6.936601 7.540327
## stream77 stream48 stream63 stream84 stream36 stream52
## 1 6.739532 6.564259 7.329269 7.255041 8.089297 9.153537
## 2 7.421621 7.321632 8.625799 8.632695 9.498303 10.266372
## 3 7.128002 7.116893 8.083832 8.221432 8.347417 7.932561
## 4 6.949546 6.892591 7.736889 7.786542 7.678440 7.309088
## 5 7.012069 6.889876 7.744266 7.698680 7.535520 7.382322toppredictions<- toppredictions[,-1]
toppredictions$upper<- 0
toppredictions$mean<- 0
toppredictions$lower<- 0
for (i in 1:nrow(toppredictions)){
# print(i)
upper= CI(as.numeric(toppredictions[i,]),ci=0.95)[1]
mean= CI(as.numeric(toppredictions[i,]),ci=0.95)[2]
lower= CI(as.numeric(toppredictions[i,]),ci=0.95)[3]
toppredictions[i,]$upper<- upper
toppredictions[i,]$mean<- mean
toppredictions[i,]$lower<- lower
}
toppredictions<-cbind(date = cwws32_run1$bd$date,toppredictions)
head(toppredictions[1:5,])## date stream58 stream19 stream126 stream125 stream57 stream137 stream40
## 1 2015-11-01 7.301086 7.621391 7.234245 7.033243 7.219342 6.935841 6.782780
## 2 2015-11-02 8.339411 8.481753 8.223325 8.076945 8.454141 8.075405 7.845375
## 3 2015-11-03 7.622130 7.759192 7.172434 7.837778 8.350942 8.026212 7.756190
## 4 2015-11-04 7.194301 7.304438 6.680168 7.257826 7.691580 7.497977 7.245324
## 5 2015-11-05 7.160050 7.311759 6.628511 7.158218 7.412052 7.315987 7.036828
## stream117 stream5 stream76 stream16 stream94 stream42 stream50 stream79
## 1 7.173283 6.999768 7.718024 7.229691 7.243884 6.920345 7.543200 6.842465
## 2 8.105602 7.945108 8.758111 8.319140 8.427707 8.165312 8.586470 7.613180
## 3 7.604787 7.850628 7.978340 7.729866 8.074425 8.206078 7.850840 7.243560
## 4 7.235594 7.204286 7.431195 7.370403 7.595946 7.660926 7.434651 6.811250
## 5 7.255528 7.144896 7.330426 7.420635 7.441549 7.419535 7.375747 6.737169
## stream128 stream80 stream127 stream146 stream51 stream49 stream8 stream62
## 1 7.031354 6.967658 6.908412 7.002993 6.982513 7.046788 7.197182 6.890182
## 2 7.979749 7.936525 7.709666 7.807816 7.912854 8.251385 8.231783 7.745624
## 3 7.685994 7.711997 7.469563 7.302992 7.773524 8.112560 7.743575 7.625534
## 4 7.303197 7.208213 6.940977 7.009744 7.280964 7.649492 7.392883 7.266316
## 5 7.224022 7.158457 6.960004 7.127518 7.109254 7.494514 7.438454 7.186921
## stream142 stream108 stream92 stream135 stream26 stream66 stream41 stream78
## 1 7.480707 7.050082 6.994545 8.149884 6.128937 7.076216 7.746484 6.794378
## 2 8.359867 7.697038 8.156216 9.069572 7.083323 8.226110 8.935218 7.660612
## 3 7.664917 7.360530 7.866016 7.740992 6.885493 7.897024 8.100632 7.610500
## 4 7.213749 7.001604 7.509209 7.264247 6.467991 7.464590 7.517883 7.036635
## 5 7.162419 7.042933 7.440129 7.260217 6.335203 7.384197 7.443158 7.061996
## stream69 stream87 stream101 stream91 stream111 stream45 stream73 stream88
## 1 7.119836 7.011195 6.720076 7.547154 6.723873 6.934614 9.046141 6.696809
## 2 7.907789 7.858544 7.854581 8.689763 7.517447 7.604283 10.119782 7.753083
## 3 7.633965 7.695322 7.900355 8.029957 6.980677 7.311125 9.198160 7.602220
## 4 7.211203 7.264788 7.409750 7.537334 6.705402 7.023489 8.598048 7.336603
## 5 7.198216 7.109192 7.301078 7.457512 6.901077 7.142527 8.370346 7.307955
## stream27 stream71 stream133 stream136 stream144 stream99 stream89 stream21
## 1 7.066230 7.179298 6.772535 9.074359 6.780157 6.669808 6.877207 6.905437
## 2 8.544112 8.333969 7.501724 9.963068 7.873735 7.499996 7.675935 8.396006
## 3 8.365502 8.008033 6.912154 7.925216 7.414847 7.176466 7.560846 8.226518
## 4 7.767062 7.515260 6.728857 7.254993 7.135348 6.913612 7.235897 7.674125
## 5 7.521184 7.458419 6.866491 7.204661 7.105699 6.960972 7.206261 7.481645
## stream28 stream107 stream67 stream22 stream81 stream110 stream122 stream20
## 1 6.645647 7.209404 6.768421 6.878164 6.894000 6.484788 7.172463 6.890570
## 2 7.596875 8.336647 7.498912 7.800314 7.756850 7.345957 8.224653 7.904469
## 3 7.715250 7.948836 7.313388 7.270891 7.206889 6.939228 7.731809 7.505268
## 4 7.161561 7.541502 7.062653 6.527299 7.006527 6.764322 7.364885 7.318275
## 5 7.145523 7.480826 7.150338 6.145675 7.169440 6.895039 7.450900 7.425878
## stream7 stream68 stream85 stream114 stream33 stream145 stream15 stream25
## 1 7.051987 7.123813 6.690226 7.013109 7.260308 6.551573 6.539392 6.454897
## 2 8.205911 8.342977 7.342435 7.974642 8.258143 7.489137 7.057671 7.339007
## 3 7.950429 8.033743 7.126135 7.487044 7.591996 7.239523 6.481352 6.869175
## 4 7.511801 7.616806 6.868980 7.244088 7.344860 7.016527 6.276343 6.729247
## 5 7.432365 7.493687 7.006537 7.407971 7.456604 7.122548 6.412394 6.899262
## stream6 stream119 stream118 stream53 stream140 stream59 stream98 stream121
## 1 6.731981 7.288448 6.697253 6.753820 6.710789 7.173417 6.622145 7.066062
## 2 7.870349 8.402888 7.489018 7.329145 7.718049 8.507835 7.299454 8.075929
## 3 8.070802 7.945091 7.234954 7.021659 7.934160 8.002741 6.976581 7.549531
## 4 7.459924 7.570510 6.919776 6.733350 7.424889 7.578688 6.747043 7.329247
## 5 7.264626 7.512899 6.909700 6.882012 7.272107 7.531734 6.941437 7.473999
## stream143 stream44 stream123 stream120 stream124 stream30 stream116 stream113
## 1 7.158570 7.370372 8.067720 6.041898 7.276809 7.248389 6.801399 6.636456
## 2 8.437932 8.183350 9.316698 6.908755 8.509898 8.522186 7.648686 7.359144
## 3 8.219739 7.664192 8.133987 6.222903 7.870076 8.105675 7.320784 7.188012
## 4 7.653680 7.728082 7.512038 6.124654 7.566975 7.698694 6.945106 6.864055
## 5 7.550356 8.447302 7.424370 6.416691 7.614377 7.597161 6.954804 6.867512
## stream109 stream139 stream102 stream65 stream47 stream83 stream134 stream14
## 1 6.606008 6.656626 7.891642 7.290588 7.798762 6.601900 7.746128 8.355276
## 2 7.289312 7.313785 9.290798 7.519508 9.089471 7.074206 8.950239 9.638516
## 3 7.038974 7.106471 8.240449 6.054104 8.014270 6.763529 8.003373 8.153606
## 4 6.837338 6.925474 7.678467 5.469140 7.481455 6.589236 7.595872 7.606671
## 5 6.936015 7.039053 7.584114 5.396264 7.433285 6.791125 7.641679 7.583419
## stream104 stream82 stream29 stream112 stream100 stream64 stream131 stream77
## 1 7.406120 7.747981 7.056639 6.655416 7.285621 6.516596 7.940796 6.739532
## 2 8.769455 8.316809 7.520984 7.387885 8.541844 7.289585 9.125171 7.421621
## 3 8.202823 7.994366 6.672738 7.050887 8.056511 7.409187 7.844993 7.128002
## 4 7.699445 7.583742 6.337906 6.710635 7.659458 6.949297 7.451631 6.949546
## 5 7.620370 7.426955 6.275115 6.741510 7.637652 6.936601 7.540327 7.012069
## stream48 stream63 stream84 stream36 stream52 upper mean lower
## 1 6.564259 7.329269 7.255041 8.089297 9.153537 7.173589 6.915854 6.658119
## 2 7.321632 8.625799 8.632695 9.498303 10.266372 8.163679 7.866137 7.568594
## 3 7.116893 8.083832 8.221432 8.347417 7.932561 7.676408 7.406455 7.136502
## 4 6.892591 7.736889 7.786542 7.678440 7.309088 7.268135 7.015501 6.762866
## 5 6.889876 7.744266 7.698680 7.535520 7.382322 7.252374 7.001687 6.751001# toppredictions[1:500,]
ggplot(toppredictions,aes(date,mean,group=1))+
geom_line(col = 'red')+
geom_ribbon(aes(ymin=lower, ymax=upper), fill = "blue", alpha=0.3)+
geom_line(data=calsubset1,aes(x=Date,y=discharge_mm))+
ylab("Discharge mm")+
xlab("Date")+
labs(title = "Top 100 Model Runs vs Measured Streamflow")
par(mfrow=c(1, 1))
{plot(paramtable$m,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - M", main = "M - Decay of hydraulic conductivity with depth")
abline(h=0,lty=2)}
{plot(paramtable$k,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - K", main = "K - Hydraulic conductivity at the surface")
abline(h=0,lty=2)}
{plot(paramtable$soil_dep,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - Sd", main = "Sd - soil depth")
abline(h=0,lty=2)}
{plot(paramtable$m_v,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - Mv", main = "Mv - Vertical decay of hydraulic conductivity with depth")
abline(h=0,lty=2)}
{plot(paramtable$k_v,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - Kv", main = "Kv - Vertical hydraulic conductivity at the surface")
abline(h=0,lty=2)}
{plot(paramtable$gw1,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - GW1", main = "GW1 - Saturated to Groundwater Coefficient")
abline(h=0,lty=2)}
{plot(paramtable$gw2,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - GW2", main = "GW2 - Groundwater Loss Coefficient")
abline(h=0,lty=2)}
old.par <- par(mfrow=c(3, 7), mar = c(5,4,4,2))
plot(paramtable$m,paramtable$NSE,xlab = "m parameter", pch = 20)
lines(lowess(paramtable$m,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$k,paramtable$NSE,xlab = "lateral ksat parameter", pch = 20)
lines(lowess(paramtable$k,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$soil_dep,paramtable$NSE,xlab = "sdepth parameter", pch = 20)
lines(lowess(paramtable$soil_dep,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$m_v,paramtable$NSE,xlab = "mv parameter", pch = 20)
lines(lowess(paramtable$m_v,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$k_v,paramtable$NSE,xlab = "ksatv parameter", pch = 20)
lines(lowess(paramtable$k_v,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$gw1,paramtable$NSE,xlab = "gw1 parameter", pch = 20)
lines(lowess(paramtable$gw1,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$gw2,paramtable$NSE,xlab = "gw2 parameter", pch = 20)
lines(lowess(paramtable$gw2,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$m,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "m parameter", pch = 20)
lines(lowess(paramtable$m,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$k,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "lateral ksat parameter", pch = 20)
lines(lowess(paramtable$k,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$soil_dep,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "sdepth parameter", pch = 20)
lines(lowess(paramtable$soil_dep,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$m_v,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "mv parameter", pch = 20)
lines(lowess(paramtable$m_v,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$k_v,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "ksatv parameter", pch = 20)
lines(lowess(paramtable$k_v,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$gw1,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "gw1 parameter", pch = 20)
lines(lowess(paramtable$gw1,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$gw2,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "gw2 parameter", pch = 20)
lines(lowess(paramtable$gw2,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$m,paramtable$KGE,ylim = c(-0.25,1),xlab = "m parameter", pch = 20)
lines(lowess(paramtable$m,paramtable$KGE, f = 0.1),col="red")
plot(paramtable$k,paramtable$KGE,ylim = c(-0.25,1),xlab = "lateral ksat parameter", pch = 20)
lines(lowess(paramtable$k,paramtable$KGE, f = 0.1),col="red")
plot(paramtable$soil_dep,paramtable$KGE,ylim = c(-0.25,1),xlab = "sdepth parameter", pch = 20)
lines(lowess(paramtable$soil_dep,paramtable$KGE, f = 0.1),col="red")
plot(paramtable$m_v,paramtable$KGE,ylim = c(-0.25,0.7),xlab = "mv parameter", pch = 20)
lines(lowess(paramtable$m_v,paramtable$KGE, f = 0.1),col="red")
plot(paramtable$k_v,paramtable$KGE,ylim = c(-0.25,0.7),xlab = "ksatv parameter", pch = 20)
lines(lowess(paramtable$k_v,paramtable$KGE, f = 0.1),col="red")
plot(paramtable$gw1,paramtable$KGE,ylim = c(-0.25,0.7),xlab = "gw1 parameter", pch = 20)
lines(lowess(paramtable$gw1,paramtable$KGE, f = 0.1),col="red")
plot(paramtable$gw2,paramtable$KGE,ylim = c(-0.25,0.7),xlab = "gw2 parameter", pch = 20)
lines(lowess(paramtable$gw2,paramtable$KGE, f = 0.1),col="red")
par(old.par)
old.par <- par(mfrow=c(1, 7), mar = c(5,4,4,2))
plot(paramtable$m,paramtable$smKGE,ylim = c(-0.25,1),xlab = "m parameter", pch = 20)
lines(lowess(paramtable$m,paramtable$smKGE, f = 0.1),col="red")
plot(paramtable$k,paramtable$smKGE,ylim = c(-0.25,1),xlab = "lateral ksat parameter", pch = 20)
lines(lowess(paramtable$k,paramtable$smKGE, f = 0.1),col="red")
plot(paramtable$soil_dep,paramtable$smKGE,ylim = c(-0.25,1),xlab = "sdepth parameter", pch = 20)
lines(lowess(paramtable$soil_dep,paramtable$smKGE, f = 0.1),col="red")
plot(paramtable$m_v,paramtable$smKGE,ylim = c(-0.25,0.7),xlab = "mv parameter", pch = 20)
lines(lowess(paramtable$m_v,paramtable$smKGE, f = 0.1),col="red")
plot(paramtable$k_v,paramtable$smKGE,ylim = c(-0.25,0.7),xlab = "ksatv parameter", pch = 20)
lines(lowess(paramtable$k_v,paramtable$smKGE, f = 0.1),col="red")
plot(paramtable$gw1,paramtable$smKGE,ylim = c(-0.25,0.7),xlab = "gw1 parameter", pch = 20)
lines(lowess(paramtable$gw1,paramtable$smKGE, f = 0.1),col="red")
plot(paramtable$gw2,paramtable$smKGE,ylim = c(-0.25,0.7),xlab = "gw2 parameter", pch = 20)
lines(lowess(paramtable$gw2,paramtable$smKGE, f = 0.1),col="red")
par(old.par)
### end plot all data for GLUE
paste0("cwws32_run",which(NSElist == Reduce(max, NSElist), arr.ind=TRUE),"$bd$streamflow")## [1] "cwws32_run142$bd$streamflow"paste0("cwws32_run",which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE),"$bd$streamflow")## [1] "cwws32_run58$bd$streamflow"paste0("cwws32_run",which(KGElist == Reduce(max, KGElist), arr.ind=TRUE),"$bd$streamflow")## [1] "cwws32_run65$bd$streamflow"{plot(cwws32_run33$bd$date,cwws32_run33$bd$streamflow, type = "l", col = 'red', main = "Observed Streamflow vs Top Runs", xlab = "Date", ylab = "Streamflow")
lines(cwws32_run65$bd$date,cwws32_run65$bd$streamflow, col = 'blue')
#lines(cwws32_run46$bd$date,cwws32_run46$bd$streamflow, col = 'green')
lines(simmerge1$Date,simmerge1$observedstreamflow, col = 'orange', lwd = 2, lty = 2)}
#{plot(paste("cwws32_run",which(NSElist == Reduce(max, NSElist), arr.ind=TRUE),"$bd$date", sep = ""),paste("cwws32_run",which(NSElist == Reduce(max, NSElist), arr.ind=TRUE),"$bd$streamflow", sep = ""), type = "l", col = 'red', main = "Observed Streamflow vs Top Runs", xlab = "Date", ylab = "Streamflow")
#lines(paste("cwws32_run",which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE),"$bd$date", sep = ""),paste("cwws32_run",which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE),"$bd$streamflow", sep = ""), col = 'blue')
#lines(cwws32_run80$bd$date,cwws32_run80$bd$streamflow, col = 'green')
#lines(paste("cwws32_run",which(KGElist == Reduce(max, KGElist), arr.ind=TRUE),"$bd$date", sep = ""),paste("cwws32_run",which(KGElist == Reduce(max, KGElist), arr.ind=TRUE),"$bd$streamflow", sep = ""), col = 'orange', lwd = 2, lty = 2)
#legend("topleft",inset = 0.05, legend=c("Oberved Discharge","Best NSE","BEST lnNSE","Best KGE"), col=c("orange","red","blue","green"), lty=c(2,1,1,1))}
{plot(cwws32_run13$bd$date,cwws32_run13$bd$streamflow, type = "l", col = 'red', main = "Observed Streamflow vs Top Runs", xlab = "Date", ylab = "Streamflow", ylim = c(0,10))
#lines(cwws32_run47$bd$date,cwws32_run47$bd$streamflow, col = 'blue')
lines(obsws32$Date,obsws32$discharge_mm, col = 'orange', lwd = 2, lty = 2)}
*** VALIDATION ***
Run Model once for validation, copied from previous code must clean up
Configure RHESSys Inputs/Outputs for a single run.
Use World Statefile created during spin-up. To run model one time at patch scale.
COMMAND LINE OPTIONS
-b Basin output option. Print out response variables for specified basins. -c Canopy stratum output option. Print out response variables for specified strata. -g Grow option. Try to read in dynamic bgc input data and output dynamic bgc parameters. -h Hillslope output option. Print out response variables for specified hillslopes. -p Patch output option. Print out response variables for specified patches. -r Routing option. Gives name of flow_table to define explicit routing connectivity. Also trigger use of explicit routing over TOPMODEL approach. -c Stratum output option. Print out response variables for specified strata.
input_rhessys = IOin_rhessys_input(
version = rh_path,
tec_file = "tecfiles/tec_daily",
world_file = "CWWS32.world.Y2018M10D31H1.state",
world_hdr_prefix = "CWWS32",
flowtable = "CWWS321.flow",
start = "2014 11 1 1",
end = "2018 11 1 1",
output_folder = "out",
output_prefix = "cwws32",
commandline_options = c("-b -g -p"))
# commandline_options = c("-b -g -p 1 128 154914 154914")
## TEC file dictates model output, begin output a year in to allow model SM to stabilize
# do not output_state or worldfile may be overwritten as output is created
input_tec_data = IOin_tec_std(start = "2015 11 1 1",
end = "2018 11 1 1",
output_state = FALSE)
input_hdr = IOin_hdr(
basin = "defs/basin.def",
hillslope = "defs/hillslope.def",
zone = "defs/zone.def",
soil = c("defs/soil_clay.def","defs/soil_clayloam.def","defs/soil_loam.def","defs/soil_loamysand.def","defs/soil_rock.def","defs/soil_sand.def","defs/soil_sandyclay.def","defs/soil_sandyclayloam.def","defs/soil_sandyloam.def","defs/soil_silt.def","defs/soil_siltyclay.def","defs/soil_siltyclayloam.def","defs/soil_siltyloam.def","defs/soil_water.def", "defs/soil_shallowloam.def", "defs/soil_shallowsandyclayloam.def", "defs/soil_shallowsandyloam.def"),
landuse = "defs/lu_undev.def",
stratum = "defs/veg_deciduous.def",
basestations = "clim/cwtws32daymet.base")
## Calibrated Values
which(NSElist == Reduce(max, NSElist), arr.ind=TRUE)## [1] 142which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE)## [1] 58which(KGElist == Reduce(max, KGElist), arr.ind=TRUE)## [1] 65which(smNSElist == Reduce(max, smNSElist), arr.ind=TRUE)## [1] 72which(smlnNSElist == Reduce(max, smlnNSElist), arr.ind=TRUE)## [1] 34which(smKGElist == Reduce(max, smKGElist), arr.ind=TRUE)## [1] 4paramtable[which(NSElist == Reduce(max, NSElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2 NSE
## 142 1.478595 6.810474 0.3167354 1.225985 14.37205 0.266005 0.1388771 0.4876977
## lnNSE KGE run smNSE smlnNSE smKGE
## 142 0.3447955 0.5987777 142 -4.221292 -2.951442 -0.2579614paramtable[which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2 NSE
## 58 1.996897 11.9516 0.1659503 14.59814 11.96051 0.04187779 0.3730755 0.4629913
## lnNSE KGE run smNSE smlnNSE smKGE
## 58 0.3570073 0.6166828 58 -5.206479 -3.536803 -0.4590018paramtable[which(KGElist == Reduce(max, KGElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2
## 65 1.669432 26.16634 0.3261013 0.3147534 11.40159 0.03675462 0.2058475
## NSE lnNSE KGE run smNSE smlnNSE smKGE
## 65 0.4219818 0.3230267 0.6270194 65 -4.489804 -3.079227 -0.1218914paramtable[which(smNSElist == Reduce(max, smNSElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2
## 72 0.8611576 1.628167 0.4270699 0.5689802 12.58125 0.2944458 0.2843299
## NSE lnNSE KGE run smNSE smlnNSE smKGE
## 72 -0.04066875 0.2492639 0.4594109 72 -1.624272 -2.020654 -0.05025428paramtable[which(smlnNSElist == Reduce(max, smlnNSElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2 NSE
## 34 1.241347 2.344012 0.4243403 13.31544 6.071543 0.1391299 0.3108362 0.09628193
## lnNSE KGE run smNSE smlnNSE smKGE
## 34 0.2993379 0.5115521 34 -1.887176 -1.50897 -0.1419453paramtable[which(smKGElist == Reduce(max, smKGElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2
## 4 0.6367556 6.600494 0.2709346 0.7792454 11.35695 0.262354 0.4609033
## NSE lnNSE KGE run smNSE smlnNSE smKGE
## 4 -0.01049478 0.1662837 0.4611981 4 -1.69059 -1.517735 0.002953977stdpars<- IOin_std_pars(
m = 0.64,
k = 6.6,
soil_dep= 0.27,
m_v = 0.78,
k_v = 11.36,
gw1 = 0.26,
gw2 = 0.46)Run RHESSys once for validation
Read RHESSys Single Run Output
getwd()## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"singlerunvalidation<- readin_rhessys_output("out/cwws32_runvalidation")
plot(singlerunvalidation$bd$streamflow~singlerunvalidation$bd$date, type = "l", main = "Streamflow", xlab = "Date", ylab = "Streamflow", col = 'DarkBlue')
plot(as.Date(singlerunvalidation$bd$date),singlerunvalidation$bd$rz_storage, type = "l", col = 'black', main = "Basin Scale Root Zone Storage",
xlab = "Date", ylab = "mm")
#basin scale soil moisture
plot(singlerunvalidation$bd$rz_storage/singlerunvalidation$bdg$root_depth~singlerunvalidation$bd$date, type = "l", main = "RZ_Storage/Root_Depth x Date", xlab = "Date", ylab = "rz_Storage/root_depth", col = 'brown')
plot(singlerunvalidation$bd$lai~singlerunvalidation$bd$date, type = "l", main = "LAI", xlab = "Date", ylab = "LAI", col = 'DarkGreen')
plot(singlerunvalidation$bd$pet~singlerunvalidation$bd$date, type = "l", main = "Potential Evapotranspiration", xlab = "Date", ylab = "PET", col = 'DarkBlue')
plot(singlerunvalidation$bd$et~singlerunvalidation$bd$date, type = "l", main = "Evapotranspiration", xlab = "Date", ylab = "ET", col = 'darkslategray')
plot((singlerunvalidation$bd$unsat_stor/singlerunvalidation$bd$sat_def_z)~singlerunvalidation$bd$date, type = "l", main = "unsat_stor/sat_def_z", xlab = "Date", ylab = "vwc", col = 'DarkGreen')
Merge Observed and Simulated Data
## merge values instead of subsetting to match up with simulated dates
singlerunvalidation$bd<- merge(singlerunvalidation$bd,obsws32, by.x = "date", by.y = "Date", all = FALSE)
## merge soil moisture values but include all data
singlerunvalidation$data <-merge(singlerunvalidation$bd,obsws32sm, by.x = "date", by.y = "Date", all = TRUE)Plot Observed and Simulated Data, Plot RHESSys Outputs
## Plot Hydrograph comparing modeled and observed streamflow
{hydrograph(input=singlerunvalidation$bd, streamflow=singlerunvalidation$bd$observedstreamflow, streamflow2 = singlerunvalidation$bd$streamflow, timeSeries = singlerunvalidation$bd$date, precip = singlerunvalidation$bd$streamflow,
P.units = "mm", S.units = "mm normalized by basin area", S1.col = 'Blue', S2.col = 'Red')
legend("topleft",inset = .2, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("blue","red"), lty=1:2)}
{plot(singlerunvalidation$data$date,singlerunvalidation$data$observedstreamflow.x, type = "l", xlab = "Date", ylab = "Streamflow", main = "Predicted vs Observed Streamflow")
lines(singlerunvalidation$data$date, singlerunvalidation$data$streamflow, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:1)}
{plot(singlerunvalidation$data$date,singlerunvalidation$data$observedstreamflow.x, type = "l", xlab = "Date", ylab = "Streamflow", main = "Predicted vs Observed Log Streamflow", log = "y")
lines(singlerunvalidation$data$date, singlerunvalidation$data$streamflow, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:1)}
predictedvsobservedstreamlm<- lm(observedstreamflow~streamflow, data = singlerunvalidation$bd)
{plot(singlerunvalidation$bd$observedstreamflow,singlerunvalidation$bd$streamflow, xlab = "Observed", ylab = "Predicted", main = "Predicted vs Observed Streamflow", xlim = c(0,55), ylim = c(0,55))
abline(a=0, b=1, col = 'RED', lty = 2, lwd = 2)
legend("right",inset = 0.01, legend = paste("R2 =",format(summary(predictedvsobservedstreamlm)$r.squared,digits=3)))}
{plot(singlerunvalidation$data$date,(singlerunvalidation$data$rz_storage/singlerunvalidation$data$rootdepth), type = "l", ylim = c(0,0.35), xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Soil Moisture")
lines(singlerunvalidation$data$date, singlerunvalidation$data$mergedsoilmoisture, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Soil Moisture","Predicted Soil Moisture"), col=c("black","red"), lty=1:1)}
predictedvsobservedsoillm<- lm((rz_storage/rootdepth)~mergedsoilmoisture, data = singlerunvalidation$data)
{plot((singlerunvalidation$data$rz_storage/singlerunvalidation$data$rootdepth),singlerunvalidation$data$mergedsoilmoisture, xlab = "Observed", ylab = "Predicted", main = "Predicted vs Observed Basin Scale Soil Moisture", xlim = c(0.05,0.4), ylim = c(0.05,0.4))
abline(a=0, b=1, col = 'RED', lty = 2, lwd = 2)
legend("right",inset = 0.01, legend = paste("R2 =",format(summary(predictedvsobservedsoillm)$r.squared,digits=3)))}
## Plot other model outputs
par(mfrow=c(5,1), mar = c(2,5,4,2))
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$tmin, type = "l", col = 'BLUE', ylim = c(-15,30), ylab = "Temperature", xlab = "Date")
lines(singlerunvalidation$bd$date,singlerunvalidation$bd$tmax, type = "l", col = 'RED')
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$precip, type = "h", ylab = "Precipitation", xlab = "Date")
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$streamflow, type = "l", ylab = "Streamflow", xlab = "Date", col = 'blue')
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$baseflow, type = "l", ylab = "Baseflow", xlab = "Date", col = 'brown')
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$rz_storage, type = "l", ylab = "RZ Storage", xlab = "Date", col = 'dark green')
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$psn ,type = "l", ylab = "PSN", xlab = "Date")
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$et, type = "l", ylab = "Evaoptranspiration", xlab = "Date")
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$pet, type = "l", ylab = "PET", xlab = "Date")
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$gw.Qout, type = "l", ylab = "Q out", xlab = "Date")
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$gw.storage, type = "l", ylab = "GW Storage", xlab = "Date")
par(mfrow=c(1,1))
### will only plot in growth mode , used to check that vegetation is initialized properly ##
plot(singlerunvalidation$bdg$date, singlerunvalidation$bdg$lai, type = "l", main = "LAI", xlab = "Date", col = 'dark green')
plot(singlerunvalidation$bdg$date, singlerunvalidation$bdg$soilc, type = "l", main = "Soil Carbon", xlab = "Date", col = 'brown')
## plot cal/val sets
{plot(singlerunvalidation$data$date,(singlerunvalidation$data$rz_storage/singlerunvalidation$data$rootdepth), type = "l", ylim = c(0,0.35), xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Soil Moisture")
lines(obsws32smcal$Date,obsws32smcal$mergedsoilmoisture, col = "BLUE")
lines(obsws32smval$Date,obsws32smval$mergedsoilmoisture, col = "RED")
legend("topright",inset = 0, legend=c("Predicted Soil Moisture","Observed Calibration Soil Moisture", "Observed Validation Soil Moisture"), col=c("black","blue","red"), lty=c(1,1,1))}
{plot(singlerunvalidation$data$date,singlerunvalidation$data$streamflow, type = "l", xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Streamflow")
lines(obsws32cal$Date,obsws32cal$observedstreamflow, col = "BLUE")
lines(obsws32val$Date,obsws32val$observedstreamflow, col = "RED")
legend("topright",inset = 0, legend=c("Predicted Streamflow","Observed Calibration Streamflow", "Observed Validation Streamflow"), col=c("black","blue","red"), lty=c(1,1,1))}
Spatial Data import and Display
# Original input maps have been cropped before processing in R, without cropping the RHESSYs preparation in the R environment will not function properly, below are the input maps used by R
ws32<- raster("grassexport/cwt_ws32/basin_ws32.tif")
patches<- raster("grassexport/cwt_ws32/patch.tif")
acc10<- raster("grassexport/cwt_ws32/acc10.tif")
aspect10<- raster("grassexport/cwt_ws32/aspect10.tif")
#bt1000<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/b.t1000.tif")
dem10<- raster("grassexport/cwt_ws32/dem10.tif")
dem10f<- raster("grassexport/cwt_ws32/dem10f.tif")
dir10<- raster("grassexport/cwt_ws32/dir10.tif")
drain10<- raster("grassexport/cwt_ws32/drain10.tif")
#ht1000<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/h.t1000.tif")
hillslope<- raster("grassexport/cwt_ws32/hillslope.tif")
roads<- brick("grassexport/cwt_ws32/roads.tif")
#shade10<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/shade10.tif")
slope10<- raster("grassexport/cwt_ws32/slope10.tif")
streams<- raster("grassexport/cwt_ws32/streams.tif")
topidx10<- raster("grassexport/cwt_ws32/topidx10_100.tif")
#xmap<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/xmap.tif")
#ymap<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/ymap.tif")
#plot(ws32, main ="ws32")
#plot(patches, main = "patches")
#plot(acc10, main = "acc10")
#plot(aspect10, main = "aspect10")
#levelplot(bt1000, col.regions = brewer.pal(name='Spectral', n = 11),at=seq(0,6500,1))
#plot(bt1000, main = "bt1000")
#plot(dem10, main = "dem10")
#plot(dem10f, main = "dem10f")
#plot(dir10, main = "dir10")
#plot(drain10, main = "drain10")
#plot(ht1000, main = "ht1000")
#plot(hillslope, main = "hillslope")
#plot(roads, main = "roads")
#plot(shade10, main = "shade10")
#plot(slope10, main = "slope10")
#plot(streams, main = "streams")
#plot(topidx10, main = "topidx10")
#plot(xmap, main = "xmap")
#plot(ymap, main = "ymap"){plot(dem10, col = colorramp, main = "Coweeta Observed Soil Moisture Locations")
points(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1]), c(ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", lwd = 1)}
matrix(c(ObservedSMLocations$Name[3],ObservedSMLocations$Name[4],ObservedSMLocations$Name[5],ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), nrow = 3, ncol = 3)## [,1] [,2] [,3]
## [1,] "132 plot center" "275432.649744231" "3881522.05364127"
## [2,] "232 plot center" "275578.214010129" "3881428.55011708"
## [3,] "332 plot center" "275760.864807773" "3881284.18476691"## Extract values for single date
displayday <- singlerunvalidation$pd[which(singlerunvalidation$pd$date=="2018-03-21")]
## reclassify all patches to na and then reclassify to model output values
reclass_all<- c(-2147483648, 2147483647,NA)
reclass_allm<- matrix(reclass_all, ncol = 3, byrow = TRUE)
displaydaypatches<- reclassify(patches,reclass_allm)
displayday$calculatedrzsm <- displayday$rz_storage/displayday$root.depth
reclass_displayday <- cbind(displayday$patchID,displayday$calculatedrzsm)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches<- reclassify(patches,reclass_displayday)
## Plot Spatial Data Output for RZ Storage
#plot(displaydaypatches, xlim = c(280900,282500), ylim = c(4870000,4872000))
{plot(mask(displaydaypatches,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 RZ Soil Moisture Prediction for 1 day", line = -2)
points(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1]), c(ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", cex = 2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}
reclass_displayday <- cbind(displayday$patchID,displayday$root.depth)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches2<- reclassify(patches,reclass_displayday)
{plot(mask(displaydaypatches2,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 root.depth Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}
reclass_displayday <- cbind(displayday$patchID,displayday$rz_storage)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches3<- reclassify(patches,reclass_displayday)
{plot(mask(displaydaypatches3,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 rz_storage Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}
reclass_displayday <- cbind(displayday$patchID,displayday$sat_def)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches4<- reclassify(patches,reclass_displayday)
{plot(mask(displaydaypatches4,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 saturation deficit Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}
## Plot Hydrograph
{hydrograph(input=singlerunvalidation$bd, streamflow=singlerunvalidation$bd$observedstreamflow, streamflow2 = singlerunvalidation$bd$streamflow, timeSeries = singlerunvalidation$bd$date, precip = singlerunvalidation$bd$streamflow,
P.units = "mm", S.units = "mm normalized by basin area", S1.col = 'Blue', S2.col = 'Red')
legend("topleft",inset = .2, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("blue","red"), lty=1:2)}
### end show patch map
patchnumbers<- raster::extract(patches,matrix(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), nrow = 3, ncol = 2))
ObserveSoilMoisturePatches <- matrix(c(ObservedSMLocations$Name[3],ObservedSMLocations$Name[4],ObservedSMLocations$Name[5],ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2],patchnumbers), nrow = 3, ncol = 4)
## Extract values for a single patch, select patch ID
ObserveSoilMoisturePatches[1,4]## [1] "136557"displaypatch<- singlerunvalidation$pd[which(singlerunvalidation$pd$patchID==ObserveSoilMoisturePatches[1,4])]
ObserveSoilMoisturePatches[2,4]## [1] "141942"displaypatch2<- singlerunvalidation$pd[which(singlerunvalidation$pd$patchID==ObserveSoilMoisturePatches[2,4])]
ObserveSoilMoisturePatches[3,4]## [1] "149478"displaypatch3<- singlerunvalidation$pd[which(singlerunvalidation$pd$patchID==ObserveSoilMoisturePatches[3,4])]
plot(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", main = 'CW WS27 More Detailed RZ Storage Patch #147322', ylim = c(0,.4))
{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", ylim = c(0,.4), col = 'darkorange', xlim=as.Date(c("2016-10-1","2019-10-1")), main = "CW WS32 More Detailed Soil Moisture vs Patches",
xlab = "Date", ylab = "Soil Moisture mm")
lines(smws32_2mean$Group.1,smws32_2mean$smois30A,type = "l", col = 'darkred')
lines(smws32_3mean$Group.1,smws32_3mean$smois30A,type = "l", col = 'blue')
lines(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'red', lty = 2)
lines(displaypatch2$date,(displaypatch2$rz_storage/displaypatch2$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'green', lty = 2)
lines(displaypatch3$date,(displaypatch3$rz_storage/displaypatch3$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'black', lty = 2)
legend("bottomleft",inset = .01, legend=c("Observed Patch 1","Predicted Patch 1","Observed Patch 2","Predicted Patch 2","Observed Patch 3","Predicted Patch 3"), col=c("darkorange", "red","darkred","green","blue","black"), lty=c(1,2,1,2,1,2))}
{plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", ylim = c(0,0.4), col = 'BLUE', xlim=as.Date(c("2016-10-1","2019-10-1")), main = "CW WS32 More Detailed Soil Moisture vs Patch #154914",
xlab = "Date", ylab = "Soil Moisture mm")
lines(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", main = 'CW WS27 More Detailed RZ Storage Patch #154914', col = 'red')
legend("bottomleft",inset = .01, legend=c("Observed Soil Moisture","Predicted Soil Moisture"), col=c("blue","red"), lty=1:1)}
#
#obsws32cal <- obsws32[obsws32$Date >= Caldates[1] & obsws32$Date <= Caldates[2], ]
#obsws32val <- obsws32[obsws32$Date >= Valdates[1] & obsws32$Date <= Valdates[2], ]
#obsws32smcal <- obsws32sm[obsws32sm$Date >= Caldates[1] & obsws32sm$Date <= Caldates[2], ]
#obsws32smval <- obsws32sm[obsws32sm$Date >= Valdates[1] & obsws32sm$Date <= Valdates[2], ]
head(obsws32cal)## Date flow discharge_mm observedstreamflow
## 22961 2017-01-01 2.30 2.172645 2.172645
## 22962 2017-01-02 3.04 2.871669 2.871669
## 22963 2017-01-03 2.72 2.569388 2.569388
## 22964 2017-01-04 2.11 1.993165 1.993165
## 22965 2017-01-05 1.92 1.813686 1.813686
## 22966 2017-01-06 1.81 1.709777 1.709777head(obsws32val)## Date flow discharge_mm observedstreamflow
## 23311 2017-12-17 2.14 2.021504 2.021504
## 23312 2017-12-18 2.27 2.144306 2.144306
## 23313 2017-12-19 2.23 2.106521 2.106521
## 23314 2017-12-20 3.39 3.202289 3.202289
## 23315 2017-12-21 2.74 2.588281 2.588281
## 23316 2017-12-22 2.60 2.456033 2.456033head(obsws32smcal)## Date flow discharge_mm observedstreamflow smois30site1a smois30site1b
## 1 2017-01-01 2.30 2.172645 2.172645 0.2194375 0.2468125
## 2 2017-01-02 3.04 2.871669 2.871669 0.2241528 0.2509896
## 3 2017-01-03 2.72 2.569388 2.569388 0.2184132 0.2409236
## 4 2017-01-04 2.11 1.993165 1.993165 0.1952743 0.2236250
## 5 2017-01-05 1.92 1.813686 1.813686 0.1825451 0.2149410
## 6 2017-01-06 1.81 1.709777 1.709777 0.1745729 0.2097500
## smois60site1a smois60site1b smois30site2a smois30site2b smois60site2a
## 1 0.1988021 0.2463715 0.2162465 0.2141181 0.2191528
## 2 0.2243646 0.2737917 0.2244410 0.2392847 0.2239028
## 3 0.2177813 0.2657361 0.2209722 0.2291771 0.2154028
## 4 0.1993090 0.2483194 0.2058437 0.2039549 0.1967431
## 5 0.1845556 0.2335799 0.1972882 0.1883021 0.1870486
## 6 0.1744132 0.2231181 0.1914861 0.1774236 0.1808646
## smois60site2b smois30site3a smois30site3b smois60site3a smois60site3b
## 1 0.1775868 0.2830729 0.2573576 0.2625312 0.2399167
## 2 0.2186528 0.2833924 0.2626007 0.2942917 0.2684410
## 3 0.2180486 0.2697604 0.2480833 0.2926944 0.2704132
## 4 0.1968403 0.2399167 0.2206285 0.2794965 0.2609514
## 5 0.1786215 0.2264861 0.2081111 0.2668299 0.2512847
## 6 0.1669340 0.2187014 0.2009896 0.2569028 0.2437292
## mergedsoilmoisture
## 1 0.2317839
## 2 0.2490255
## 3 0.2422839
## 4 0.2225752
## 5 0.2099661
## 6 0.2015738head(obsws32smval)## Date flow discharge_mm observedstreamflow smois30site1a smois30site1b
## 351 2017-12-17 2.14 2.021504 2.021504 0.1181493 0.2395313
## 352 2017-12-18 2.27 2.144306 2.144306 0.1165625 0.2392708
## 353 2017-12-19 2.23 2.106521 2.106521 0.1135139 0.2332326
## 354 2017-12-20 3.39 3.202289 3.202289 0.1484514 0.2609861
## 355 2017-12-21 2.74 2.588281 2.588281 0.1391528 0.2532083
## 356 2017-12-22 2.60 2.456033 2.456033 0.1312049 0.2428194
## smois60site1a smois60site1b smois30site2a smois30site2b smois60site2a
## 351 0.1677882 0.2054826 0.1575382 0.1352778 NA
## 352 0.1635625 0.2064132 0.1605903 0.1370000 NA
## 353 0.1600174 0.2057812 0.1585660 0.1363750 NA
## 354 0.1844479 0.2319132 0.1872778 0.1721458 NA
## 355 0.2092639 0.2517917 0.1925139 0.1736111 NA
## 356 0.1954965 0.2400937 0.1847604 0.1649479 NA
## smois60site2b smois30site3a smois30site3b smois60site3a smois60site3b
## 351 0.09663194 0.2347188 0.1869757 0.2350000 0.2169896
## 352 0.09784375 0.2385694 0.1832847 0.2350000 0.2167847
## 353 0.09828472 0.2289757 0.1792153 0.2350000 0.2160000
## 354 0.12405208 0.2609826 0.2203681 0.2421354 0.2173993
## 355 0.14059028 0.2525278 0.2136979 0.2712326 0.2365937
## 356 0.13011806 0.2397049 0.1995382 0.2683160 0.2404201
## mergedsoilmoisture
## 351 0.1812803
## 352 0.1813529
## 353 0.1786329
## 354 0.2045600
## 355 0.2121985
## 356 0.2034018#subset Validation range based on Val Dates set by available SM data
validationsubset <- singlerunvalidation$bd[singlerunvalidation$bd$date >= Valdates[1] & singlerunvalidation$bd$date <= Valdates[2], ]
## Fit tests for streamflow
NSE(validationsubset$streamflow,validationsubset$observedstreamflow)## [1] 0.151659NSE(validationsubset$streamflow,validationsubset$observedstreamflow, FUN = log, epsilon = "Pushpalatha2012", na.rm=TRUE)## [1] 0.6005506KGE(validationsubset$streamflow,validationsubset$observedstreamflow)## [1] 0.5055394# Plot Observed and Predicted Streamflow
{plot(validationsubset$date,validationsubset$streamflow, type = 'l', lty = 2, col = 'RED', main = "Validation Subset Observed and Predicted Streamflow")
lines(validationsubset$date,validationsubset$observedstreamflow, type = 'l', col = 'BLACK')
legend("topright",inset = .1, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:2)}
{plot(singlerunvalidation$bd$date,singlerunvalidation$bd$streamflow, type = 'l', lty = 2, col = 'RED', main = "Observed and Predicted Streamflow")
lines(singlerunvalidation$bd$date,singlerunvalidation$bd$observedstreamflow, type = 'l', col = 'BLACK')
legend("topright",inset = .1, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:2)}
#Plot Model RZ SM outputs for patch
{plot(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l",ylim = c(0,0.50), ylab = "Root Zone Storage (mm/m) ", xlab = "Date", main = 'Root Zone Display Patch')
lines(displaypatch$date,(displaypatch$rz_field_capacity/displaypatch$root.depth), type = "l", col = "blue")
lines(displaypatch$date,(displaypatch$rz_wilting_point/displaypatch$root.depth), type = "l", col = "red")
legend("bottomleft",inset = .1, legend=c("Predicted RZ Storage","RZ Field Capacity","RZ Wilting Point"), col=c("black","blue","red"), lty=1:1)}
Model Validation
Valdates## [1] "2017-12-17" "2018-05-15"## for now run top run and record values
paramtable[which(KGElist == Reduce(max, KGElist), arr.ind=TRUE),]## m k soil_dep m_v k_v gw1 gw2
## 65 1.669432 26.16634 0.3261013 0.3147534 11.40159 0.03675462 0.2058475
## NSE lnNSE KGE run smNSE smlnNSE smKGE
## 65 0.4219818 0.3230267 0.6270194 65 -4.489804 -3.079227 -0.1218914{plot(cwws32_run4$bd$date,cwws32_run4$bd$streamflow, type = "l", col = 'red', main = "Observed Streamflow vs Top KGE Run", xlab = "Date", ylab = "Streamflow")
lines(obsws32$Date,obsws32$discharge_mm, col = 'orange', lwd = 2, lty = 2)}
valsubset <- simmerge4[simmerge4$Date >= Valdates[1] & simmerge4$Date <= Valdates[2], ]
ValKGE<- KGE(sim = valsubset$streamflow, obs = valsubset$observedstreamflow)
valsubsetsm<- simmergesm4[simmergesm4$Date >= Valdates[1] & simmergesm4$Date <= Valdates[2],]
ValsmKGE<- KGE(sim = valsubsetsm$rz_storage/valsubsetsm$rootdepth, obs = valsubsetsm$mergedsoilmoisture)
ValKGE## [1] 0.5061636ValsmKGE## [1] -0.610026